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The emergence of synchronization in a network of coupled oscilla- 
tors is a fascinating topic in various scientific disciplines. A coupled 
oscillator network is characterized by a population of heterogeneous 
oscillators and a graph describing the interaction among them. It is 
known that a strongly coupled and sufficiently homogeneous network 
synchronizes, but the exact threshold from incoherence to synchrony 
is unknown. Here we present a novel, concise, and closed-form con- 
dition for synchronization of the fully nonlinear, non-equilibrium, and 
dynamic network. Our synchronization condition can be stated ele- 
gantly in terms of the network topology and parameters, or equiv- 
alently in terms of an intuitive, linear, and static auxiliary system. 
Our results significantly improve upon the existing conditions advo- 
cated thus far, they are provably exact for various interesting network 
topologies and parameters, they are statistically correct for almost 
all networks, and they can be applied equally to synchronization 
phenomena arising in physics and biology as well as in engineered 
oscillator networks such as electric power networks. We illustrate the 
validity, the accuracy, and the practical applicability of our results in 
complex networks scenarios and in smart grid applications. 
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The scientific interest in the synchronization of coupled 
oscillators can be traced back to Christiaan Huygens' seminal 
work on "an odd kind sympathy" between coupled pendulum 
clocks [1], and it continues to fascinate the scientific commu- 
nity to date [2, 3]. A mechanical analog of a coupled oscillator 
network is shown in Figure 1 and consists of a group of parti- 
cles constrained to rotate around a circle and assumed to move 
without colliding. Each particle is characterized by a phase 
angle 0i and has a preferred natural rotation frequency uji. 
Pairs of interacting particles i and j are coupled through an 
elastic spring with stiffness cnj . Intuitively, a weakly coupled 
oscillator network with strongly heterogeneous natural fre- 
quencies uji does not display any coherent behavior, whereas a 
strongly coupled network with sufficiently homogeneous nat- 
ural frequencies is amenable to synchronization. These two 
qualitatively distinct regimes are illustrated in Figure 1. 




Fig. 1. Mechanical analog of a coupled oscillator network (a) and its dynamics in a 
strongly coupled (b) and weakly coupled (c) network. With exception of the coupling 
weights dij, all parameters in the simulation (b) and (c) are identical. 



Formally, the interaction among n such phase oscilla- 
tors is modeled by a connected graph G(V,£,A) with nodes 



V = {l,...,n}, edges S C V x V, and positive weights 
dij > for each undirected edge £ S. For pairs of non- 

interacting oscillators i and j, the coupling weight aij is zero. 
We assume that the node set is partitioned as V = Vi U V2, and 
we consider the following general coupled oscillator model: 

Midi + DA =Ui- V n _ (Hj sin(0< - 0j) , ieVi, 

[1] 

DiOi = uji - ._ l cnj sin(0i — 6j) , i e V2 . 

The coupled oscillator model [ 1 ] consists of the second-order 
oscillators Vi with Newtonian dynamics, inertia coefficients 
Mi, and viscous damping Di. The remaining oscillators V2 
feature first-order dynamics with time constants Di . A perfect 
electrical analog of the coupled oscillator model [ 1 ] is given by 
the classic structure-preserving power network model [4], our 
enabling application of interest. Here, the first and second- 
order dynamics correspond to loads and generators, respec- 
tively, and the right-hand sides depict the power injections uji 
and the power flows aij sin(#i — 0j) along transmission lines. 

The rich dynamic behavior of the coupled oscillator 
model [1] arises from a competition between each oscilla- 
tor's tendency to align with its natural frequency uji and 
the synchronization-enforcing coupling aij sin(0i —0j) with its 
neighbors. In absence of the first term, the coupled oscillator 
dynamics [ 1 ] collapse to a trivial phase-synchronized equilib- 
rium, where all angles 0i are aligned. The dissimilar natural 
frequencies Ui, on the other hand, drive the oscillator network 
away from this all-aligned equilibrium. Moreover, even if the 
coupled oscillator model [1] synchronizes, it still carries the 
flux of angular rotation, respectively, the flux of electric power 
from generators to loads in a power network. The main and 
somehow surprising result of this paper is that, in spite of 
all the aforementioned complications, an elegant and easy to 
verify criterion characterizes synchronization of the nonlinear 
and non-equilibrium dynamic oscillator network [1]. 

Review of Synchronization in Oscillator Networks 

The coupled oscillator model [1] unifies various models in 
the literature including dynamic models of electric power net- 
works. The supplementary information (SI) discusses mod- 
eling of electric power networks in detail. For V2 = 0, the 
coupled oscillator model [1] appears in synchronization phe- 
nomena in animal flocking behavior [5] , populations of flashing 
fireflies [6], crowd synchrony on London's Millennium bridge 
[7], as well as Huygen's pendulum clocks [8]. For Vi = 0, the 
coupled oscillator model [1] reduces to the celebrated Ku- 
ramoto model [9], which appears in coupled Josephson junc- 
tions [10], particle coordination [11], spin glass models [12, 13], 
neuroscience [14], deep brain stimulation [15], chemical oscil- 
lations [16], biological locomotion [17], rhythmic applause [18], 
and countless other synchronization phenomena [19, 20, 21]. 
Finally, coupled oscillator models of the form [ 1 ] also serve as 
prototypical examples in complex networks studies [22, 23]. 

The coupled oscillator dynamics [1] feature the synchro- 
nizing effect of the coupling described by the graph G(V, £, A) 



and the de-synchronizing effect of the dissimilar natural fre- 
quencies LJi. The complex network community asks questions 
of the form "what are the conditions on the coupling and the 
dissimilarity such that a synchronizing behavior emerges?" 
Similar questions appear also in all the aforementioned ap- 
plications, for instance, in large-scale electric power systems. 
Since synchronization is pervasive in the operation of an in- 
terconnected power grid, a central question is "under which 
conditions on the network parameters and topology, the cur- 
rent load profile and power generation, does there exist a syn- 
chronous operating point [24, 25], when is it optimal [26], 
when is it stable [27, 28], and how robust is it [29, 37, 31, 32]?" 
A local loss of synchrony can trigger cascading failures and 
possibly result in wide-spread blackouts. In the face of the 
complexity of future smart grids and the integration chal- 
lenges posed by renewable energy sources, a deeper under- 
standing of synchronization is increasingly important. 

Despite the vast scientific interest, the search for sharp, 
concise, and closed- form synchronization conditions for cou- 
pled oscillator models of the form [1] has been so far in vain. 
Loosely speaking, synchronization occurs when the coupling 
dominates the dissimilarity. Various conditions have been pro- 
posed to quantify this trade-off [21, 32, 28, 33, 22, 23, 31, 34]. 
The coupling is typically quantified by the nodal degree or the 
algebraic connectivity of the graph G, and the dissimilarity 
is quantified by the magnitude or the spread of the natural 
frequencies uji. Sometimes, these conditions can be evalu- 
ated only numerically since they depend on the network state 
[32, 31] or arise from a non-trivial linearization process, such 
as the Master stability function formalism [22, 23]. To date, 
exact synchronization conditions are known only for simple 
coupling topologies [17, 21, 35, 36]. For arbitrary topologies 
only sufficient conditions are known [32, 28, 33, 31] as well 
as numerical investigations for random networks [37, 38, 39]. 
Simulation studies indicate that the known sufficient condi- 
tions are very conservative estimates on the threshold from 
incoherence to synchrony. Literally, every review article on 
synchronization concludes emphasizing the quest for exact 
synchronization conditions for arbitrary network topologies 
and parameters [20, 21, 19, 22, 23]. In this article, we present 
a concise and sharp synchronization condition which features 
elegant graph-theoretic and physical interpretations. 



Novel Synchronization Condition 

For the coupled oscillator model [1] and its applications, the 
following notions of synchronization are appropriate. First, a 
solution has synchronized frequencies if all frequencies Si are 
identical to a common constant value cj syn c. If a synchronized 
solution exists, it is known that the synchronization frequency 
is cJsync = X}fc=i u k/ D k and that, by working in a ro- 
tating reference frame, one may assume cj sy nc = 0. Second, a 
solution has cohesive phases if every pair of connected oscilla- 
tors has phase distance smaller than some angle 7 £ [0, 7r/2[, 
that is, \0i — 0j\ < 7 for every edge {i,j} G S. 

Based on a novel analysis approach to the synchronization 
problem, we propose the following synchronization condition 
for the coupled oscillator model [ 1 ] : 

Sync condition: The coupled oscillator model [1] has 
a unique and stable solution #* with synchronized fre- 
quencies and cohesive phases \0* — 9* \ < 7 < tt/2 for 
every pair of connected oscillators {i,j} G £ if 

||L^|| £oo <sin(7). [2] 



Here, Lr is the pseudo-inverse of the network Laplacian 
matrix L and \\x\\ £ ^ — max^}^ \xi—Xj \ is the worst- 
case dissimilarity for x = (xi, . . . , x n ) over the edges S. 

We establish the broad applicability of the proposed condition 
[2] to various classes of networks via analytical and statisti- 
cal methods in the next section. Before that, we provide some 
equivalent formulations for condition [2] in order to develop 
deeper intuition and obtain insightful conclusions. 

Complex network interpretation: Surprisingly, topo- 
logical or spectral connectivity measures such as nodal degree 
or algebraic connectivity are not key to synchronization. In 
fact, these often advocated [32, 28, 33, 31, 22, 23] connec- 
tivity measures turn out to be conservative estimates of the 
synchronization condition [2]. This statement can be seen 
by introducing the matrix U of orthonormal eigenvectors of 
the network Laplacian matrix L with corresponding eigenval- 
ues = Ai < A2 < • • • < A n . From this spectral viewpoint, 
condition [2] can be equivalently written as 

||t/diag(0,l/A 2 ,...,l/A n ) • (U T w)\\ e >oo <sin( 7 ). [3] 

In words, the natural frequencies uj are projected on the net- 
work modes U, weighted by the inverse Laplacian eigenval- 
ues, and || • \\s,oo evaluates the worst-case dissimilarity of 
this weighted projection. A sufficient condition for the in- 
equality [3] to be true is the algebraic connectivity condition 
A2 > I Ml £,00 • sin(7). Likewise, a necessary condition for 
inequality [3] is 2 • deg(C) > A n > ||cj||£,oo • sin (7), where 
deg(G) is the maximum nodal degree in the graph G(V, £, A). 
Clearly, when compared to [3], this sufficient condition and 
this necessary condition feature only one of n — 1 non-zero 
Laplacian eigenvalues and are overly conservative. 

Kuramoto oscillator perspective: Notice, that in the 
limit 7 — > 7r/2, condition [2] suggests that there exists a sta- 
ble synchronized solution if 

Wu\\ P <1. [4] 
11 lit, 00 L J 

For classic Kuramoto oscillators coupled in a complete graph 
with uniform weights = K/n, the synchronization condi- 
tion [4] reduces to the condition K > maXj jG { lr--;n } \oJi~ ujj\, 
known for the classic Kuramoto model [21]. 

Power network perspective: In power systems engi- 
neering, the equilibrium equations of the coupled oscillator 
model [1], given by uji — Sj=i a ij sm (^ are referred to 

as the AC power flow equations, and they are often approxi- 
mated by their linearization [29, 30, 31, 32] uji — Sj=i a ij(^i~ 
0j), known as the DC power flow equations. In vector notation 
the DC power flow equations read as uj — L6, and their solu- 
tion satisfies max{^} e £ \0i — 0j\ = \\L^u\\£ )0 o- According to 
condition [2], the worst phase distance ||L^cj||£,oo obtained by 
the DC power flow equations needs to be smaller than sin(7), 
such that the solution to the AC power flow equations satis- 
fies max{ i; j} G £: \0i —0j \ < 7. Hence, our condition extends the 
common DC power flow approximation from infinitesimally 
small angles 7 <C 1 to large angles 7 £ [0, 7r/2[. 

Auxiliary linear perspective: As detailed in the previ- 
ous paragraph, the key term L^uj in condition [2] equals the 
phase differences obtained by the linear Laplacian equation 
uj = LB. This linear interpretation is not only insightful but 
also practical since condition [2] can be quickly evaluated by 
numerically solving the sparse linear system uj — L0. Despite 
this linear interpretation, we emphasize that our derivation of 
condition [2] is not based on any linearization arguments. 

Energy landscape perspective: Condition [2] can also 
be understood in terms of an appealing energy landscape in- 
terpretation. The coupled oscillator model [1] is a system of 
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particles that aim to minimize the energy function 

£ w = E {iJ}6£ i 1 - cos ^ - 9 i)) - EL Wi ■ 6i > 

where the first term is a pair-wise nonlinear attraction among 
the particles, and the second term represents the external 
force driving the particles away from the "all-aligned" state. 
Since the energy function E(6) is difficult to study, it is nat- 
ural to look for a minimum of its second-order approximation 
Eo(0) = Ys{i,j}e£ a ij(Qi - #j) 2 / 2 ~ ELi^ ' where the 
first term corresponds to a Hookean potential. Condition [2] 
is then restated as follows: E(0) features a phase cohesive 
minimum with interacting particles no further than 7 apart if 
Eq(6) features a minimum with interacting particles no fur- 
ther from each other than sin(7), as illustrated in Figure 2. 




Fig. 2. The energy function E(0) and its quadratic approximation Eq{6) for a 
two-particle system are shown as solid and dashed curves, respectively, for the stable 
(blue), marginal (green) and unstable (red) cases. The circles and diamonds represent 
stable critical points of E(9) and Eq(0). 



Statistical result: With 99.97 % probability, for a 
nominal network, condition [2] guarantees the existence 
of an unique and stable solution #* with synchronized 
frequencies and cohesive phases \6* — 0j \ < 7 for every 
pair of connected oscillators {i, j} £ S. 

From this statistical result, we deduce that the proposed syn- 
chronization condition [2] holds true for almost all network 
topologies and parameters. Indeed, we also show the existence 
of possibly-thin sets of topologies and parameters for which 
our condition [2] is not sufficiently tight. We refer to the SI 
for an explicit family of carefully engineered and "degener- 
ate" counterexamples. Overall, our analytical and statistical 
results validate the correctness of the proposed condition [4]. 

After having established the statistical correctness of con- 
dition [2], we now investigate its predictive power for arbi- 
trary networks. Since we analytically establish that condi- 
tion [2] is exact for sufficiently small pairwise phase cohe- 
siveness \0i — 0j\ <C 1, we now investigate the other extreme, 
max/ i; j} G £: \6i — 0j\ = tt/2. To test the corresponding condi- 
tion [4] in a low-dimensional parameter space, we consider a 
complex network of Kuramoto oscillators 

6i=uJi-K-y dij sin(0i - 6j) , i G {1, . . . , n} , [5] 
z — 'j=i 

where all coupling weights dij are either zero or one, and the 
coupling gain K > serves as control parameter. If L is the 
corresponding unweighted Laplacian matrix, then condition 
[4] reads as K > ^critical = \\L (jJ £.00 • Of course, the con- 
dition K > identical is only sufficient and the critical coupling 



Analytical and Statistical Results 

Our analysis approach to the synchronization problem is 
based on algebraic graph theory. We propose an equivalent re- 
formulation of the synchronization problem, which reveals the 
crucial role of cycles and cut-sets in the graph and ultimately 
leads to the synchronization condition [2]. In particular, we 
analytically establish the synchronization condition [2] for 
the following six interesting cases: 

Analytical result: The synchronization condition [2] 
is necessary and sufficient for (i) the sparsest (acyclic) 
and (ii) the densest (complete and uniformly weighted) 
network topologies C(V, S, A), (Hi) the best (phase syn- 
chronizing) and (iv) the worst (cut-set inducing) natu- 
ral frequencies, (v) for cyclic topologies of length strictly 
less than five, (vi) for arbitrary cycles with symmetric 
parameters, (vii) as well as one-connected combinations 
of networks each satisfying one of the conditions (i)-(vi). 

A detailed and rigorous mathematical derivation and state- 
ment of the above analytical result can be found in the SI. 

After having analytically established condition [2] for a 
variety of particular network topologies and parameters, we 
establish its correctness and predictive power for arbitrary 
networks. Extensive simulation studies lead to the conclusion 
that the proposed synchronization condition [2] is statisti- 
cally correct. In order to verify this hypothesis, we conducted 
Monte Carlo simulation studies over a wide range of natural 
frequencies uji, network sizes n, coupling weights dij, and dif- 
ferent random graph models of varying degrees of sparsity and 
randomness. In total, we constructed 1.2- 10 6 samples of nomi- 
nal random networks, each with a connected graph G(V, £, A) 
and natural frequencies uj satisfying HL^cjU^oo < sin (7) for 
some 7 < 7r/2. The detailed results can be found in the SI 
and allow us to establish the following probabilistic result with 
a confidence level of at least 99% and accuracy of at least 99%: 
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Fig. 3. Numerical evaluation of the exact critical coupling K in a complex Ku- 
ramoto oscillator network. The subfigures show K normalized by ||L^o;||£ )00 for an 
Erdos-Renyi graph with probability p of connecting two nodes, for a random geometric 
graph with connectivity radius p, and for a Watts-Strogatz small world network with 
rewiring probability p. Each data point is the mean over 100 samples of the respective 
random graph model, for values of OJi sampled from a bipolar or a uniform distribution 
supported on [—1, 1], and for the network sizes n G {10, 20, 40, 80, 160}. 
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may be smaller than if critical • In order to test the accuracy of 
the condition K > if critical , we numerically found the smallest 
value of K leading to synchrony with phase cohesiveness tt/2. 

Figure reports our findings for various network sizes, con- 
nected random graph models, and sample distributions of the 
natural frequencies. We refer to the SI for the detailed simu- 
lation setup. First, notice from Subngures (a),(b),(d), and (e) 
that condition [4] is extremely accurate for a sparse graph, 
that is, for small p and n, as expected from our analyti- 
cal results. Second, for a dense graph with p « 1, Subng- 
ures (a),(b),(d), and (e) confirm the results known for clas- 
sic Kuramoto oscillators [21]: for a bipolar distribution con- 
dition [4] is exact, and for a uniform distribution a small 
critical coupling is obtained. Third, Subngures (c) and (d) 
show that condition [4] is scale- free for a Watts-Strogatz 
small world network, that is, it has almost constant accu- 
racy for various values of n and p. Fourth and finally, ob- 
serve that condition [4] is always within a constant factor of 
the exact critical coupling, whereas other proposed conditions 
[32, 28, 33, 31, 22, 23] on the nodal degree or on the algebraic 
connectivity scale poorly with respect to network size n. 



Applications in Power Networks 

We envision that condition [2] can be applied to quickly as- 
sess synchronization and robustness in power networks un- 
der volatile operating conditions. Since real- world power net- 
works are carefully engineered systems with particular net- 
work topologies and parameters, we do not extrapolate the 
statistical results from the previous section to power grids. 
Rather, we consider ten widely-established IEEE power net- 
work test cases provided by [40, 41]. 

Under nominal operating conditions, the power generation 
is optimized to meet the forecast demand, while obeying the 
AC power flow laws and respecting the thermal limits of each 
transmission line. Thermal limits constraints are precisely 
equivalent to phase cohesiveness requirements. In order to 
test the synchronization condition [2] in a volatile smart grid 
scenario, we make the following changes to the nominal net- 
work: 1) We assume fluctuating demand and randomize 50% 
of all loads to deviate from the forecasted loads. 2) We assume 
that the grid is penetrated by renewables with severely fluc- 
tuating power outputs, for example, wind or solar farms, and 
we randomize 33% of all generating units to deviate from the 
nominally scheduled generation. 3) Following the paradigm 




of smart operation of smart grids [42] , the fluctuations can be 
mitigated by fast-ramping generation, such as fast-response 
energy storage including batteries and flywheels, and control- 
lable loads, such as large-scale server farms or fleets of plug-in 
hybrid electrical vehicles. Here, we assume that the grid is 
equipped with 10% fast-ramping generation and 10% control- 
lable loads, and the power imbalance (caused by fluctuating 
demand and generation) is uniformly dispatched among these 
adjustable power sources. For each of the ten IEEE test cases, 
we construct 1000 random realizations of the scenario 1), 2), 
and 3) described above, we numerically check for the exis- 
tence of a synchronous solution, and we compare the numeri- 
cal solution with the results predicted by our synchronization 
condition [2]. Our findings are reported in Table 2, and a 
detailed description of the simulation setup can be found in 
the SI. It can be observed that condition [2] predicts the cor- 
rect phase cohesiveness \0i — 0j\ along all transmission lines 
G S with extremely high accuracy even for large-scale 
networks featuring 2383 nodes. 

As a final test, we validate the synchronization condi- 
tion [2] in a stressed power grid case study. We consider 
the IEEE Reliability Test System 96 (RTS 96) [41] illustrated 
in Figure . We assume the following two contingencies have 
taken place and we characterize the remaining safety mar- 
gin. First, we assume generator 323 is disconnected, possibly 
due to maintenance or failure events. Second, we consider 
the following imbalanced power dispatch situation: the power 
demand at each load in the Southeastern area deviates from 
the nominally forecasted demand by a uniform and positive 
amount, and the resulting power deficiency is compensated 
by uniformly increasing the generation in the Northwestern 
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Fig. 4. Illustration of contingencies the RTS 96 power network. Here, square 
nodes are generators and round nodes are loads, large amounts of power are exported 
from the Northwestern area to the Southeastern area, and generator 323 is tripped. 
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Fig. 5. The RTS 96 dynamics for a continuous load increase from 22.19% to 
22.24%. Subfigure (a) shows the angles 6(t) which loose synchrony at t* = 
18.94 s, when the thermal limit 7* = 0.1977 rad of the transmission line 
{121,325} is reached. Subfigure (b) shows the angles 9(t) at t = t* . Sub- 
figure (c) depicts the angular distances and the thermal limits 7* and 7** of the 
lines {121, 325} and {223, 318}. Subfigures (d) and (e) show the generator phase 
space (0(t),0(t)) before and after t* , where the loss of a common synchronization 
frequency can be observed. 
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area. This imbalance can arise, for example, due to a short- 
fall in predicted load and renewable energy generation. Cor- 
respondingly, power is exported from the Northwestern to the 
Southeastern area via the transmission lines {121,325} and 
{223,318}. At a nominal operating condition, the RTS 96 
power network is sufficiently robust to tolerate each single 
one of these two contingencies, but the safety margin is now 
minimal. When both contingencies are combined, then our 
synchronization condition [2] predicts that the thermal limit 
of the transmission line {121, 325} is reached at an additional 
loading of 22.20%. Indeed, the dynamic simulation scenario 
shown in Figure validates the accuracy of this prediction. It 
can be observed, that synchronization is lost for an additional 
loading of 22.33%, and the areas separate via the transmission 
line {121,325}. This separation triggers a cascade of events, 
such as the outage of the transmission line {223, 318}, and the 
power network is en route to a blackout. We remark that, if 
generator 323 is not disconnected and there are no thermal 
limit constraints, then, by increasing the loading, we observe 
the classic loss of synchrony through a saddle-node bifurca- 
tion. Also this bifurcation can be predicted accurately by our 
results, see the SI for a detailed description. 

In summary, the results in this section confirm the valid- 
ity, the applicability, and the accuracy of the synchronization 
condition [2] in complex power network scenarios. 



Discussion and Conclusions 

In this article we studied the synchronization phenomenon 
for broad class of coupled oscillator models proposed in the 
scientific literature. We proposed a surprisingly simple condi- 
tion that accurately predicts synchronization as a function of 
the parameters and the topology of the underlying network. 
Our result, with its physical and graph theoretical interpre- 
tations, significantly improves upon the existing test in the 
literature on synchronization. The correctness of our syn- 
chronization condition is established analytically for various 
interesting network topologies and via Monte Carlo simula- 
tions for a broad range of generic networks. We validated our 
theoretical results for complex Kuramoto oscillator networks 
as well as in smart grid applications. 

Our results equally answer as many questions as they pose. 
Among the important theoretical problems to be addressed is 
a characterization of the set of all network topologies and pa- 
rameters for which our proposed synchronization condition 
||^^||^,oo < 1 is not sufficiently tight. We conjecture that 



this set is "thin" in an appropriate parameter space. Our 
results suggest that an exact condition for synchronization 
of any arbitrary network is of the form ||L^o;||£ j00 < c, and 
we conjecture that the constant c is always strictly positive, 
upper-bounded, and close to one. Yet another important ques- 
tion not addressed in the present article concerns the region 
of attraction of a synchronized solution. We conjecture that 
the latter depends on the gap in the presented synchroniza- 
tion condition. On the application side, we envision that our 
synchronization conditions enable emerging smart grid appli- 
cations, such as power flow optimization subject to stability 
constraints, distance to failure metric, and the design of con- 
trol strategies to avoid cascading failures. 



Table 1. Evaluation of condition [2] for ten IEEE test 
cases under volatile operating conditions. 



Randomized test 
case (1000 instances): 


* Correctness: 


^Accuracy: 


* Cohesive 
phases: 


Chow 9 bus system 


always true 


4.1218 • lO" 5 


0.12889 


IEEE 14 bus system 


always true 


2.7995 • lO" 4 


0.16622 


IEEE RTS 24 


always true 


1.7089 • 1Q-' 6 


0.22309 


IEEE 30 bus system 


always true 


2.6140 • lO" 4 


0.1643 


New England 39 


always true 


6.6355 • lO" 5 


0.16821 


IEEE 57 bus system 


always true 


2.0630 • 10~ 2 


0.20295 


IEEE RTS 96 


always true 


2.6076 • 10~ a 


0.24593 


IEEE 118 bus system 


always true 


5.9959 • lO" 4 


0.23524 


IEEE 300 bus system 


always true 


5.2618 • lO" 4 


0.43204 


Polish 2383 bus 
system (winter 99) 


always true 


4.2183 • 10~ a 


0.25144 



Correctness: ||L+o;||£ >00 < sin(7) max^ ^g^ |0* - 0* | < 7 

t Accuracy: max^jj^ \0* - * | - arcsin(||Lt^|| £)00 ) 
* Phase cohesiveness: max{^j e £ \6* — 9 j \ 



^The accuracy and phase cohesiveness results in the third and 
fourth column are given in the unit [rad], and they are averaged 
over 1000 instances of randomized load and generation. 
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Supplementary Information 



Introduction 

This supplementary information is organized as follows. 

The section Mathematical Models and Synchronization 
Notions provides a description of the considered coupled os- 
cillator model including a detailed modeling of a mechanical 
analog and a few power network models. Furthermore, we 
state our definition of synchronization and compare various 
synchronization conditions proposed for oscillator networks. 

The section Mathematical Analysis of Synchronization 
provides a rigorous mathematical analysis of synchronization, 
which leads to the novel synchronization conditions proposed 
in the main article. Throughout our analysis we provide vari- 
ous examples illustrating certain theoretical concepts and re- 
sults, and we also compare our results to existing results in 
the synchronization and power networks literature. 

The section Statistical Synchronization Assessment pro- 
vides a detailed account of our Monte Carlo simulation studies 
and the complex Kuramoto network studies. Throughout this 
section, we also recall the basics of probability estimation by 
Monte Carlo methods that allow us to establish a statistical 
synchronization result in a mathematically rigorous way. 

Finally, the section Synchronization Assessment for Power 
Networks describes the detailed simulation setup for the ran- 
domized IEEE test systems, it provides the simulation data 
used for the dynamic IEEE RTS 96 power network simula- 
tions, and it illustrates a dynamic bifurcation scenario in the 
IEEE RTS 96 power network. 

The remainder of this section introduces some notation 
and recalls some preliminaries. 

Preliminaries and Notation. Vectors and functions: Let l n 
and n be the n-dimensional vector of unit and zero en- 
tries, and let 1^ be the orthogonal complement of l n in R n , 
that is, 1^; = {x £ R n : x _L l n }. Let be zth canoni- 
cal basis vector of R n , that is, the zth entry of ef is 1 and 
all other entries are zero. Given an n-tuple (xi,...,x n ), 
let x £ R n be the associated vector. For an ordered in- 
dex set X of cardinality \L\ and an one-dimensional array 
we define diag({c;}; e z) £ R' :r ' x ' :z: ' to be the asso- 
ciated diagonal matrix. For x £ R n , define the vector- valued 
functions sin(x) = (sin(xi), . . . , sin(x n )) and arcsin(x) = 
(arcsin(xi), . . . , arcsin(x n )), where the arcsin function is de- 
fined for the branch [-7r/2, tt/2]. For a set X C R n and a 
matrix A £ R mXn , let AX = {y £ R m : y = Ax , x £ X}. 

Geometry on n-torus: The set S 1 denotes the unit circle, 
an angle is a point £ S 1 , and an arc is a connected subset 
of S 1 . The geodesic distance between two angles 0\, O2 £ S 1 
is the minimum of the counter-clockwise and the clockwise 
arc length connecting 0\ and 62- With slight abuse of nota- 
tion, let 1 0i — O2 1 denote the geodesic distance between two 
angles 61,62 £ S 1 . Finally, the n-torus is the product set 
T n = S 1 x • • • x S 1 is the direct sum of n unit circles. 

Algebraic graph theory: Given an undirected, connected, 
and weighted graph G(V, £, A) induced by the symmet- 
ric, irreducible, and nonnegative adjacency matrix A £ 
R nXn , the Laplacian matrix L £ R nXn is defined by L = 
dia g({E^i a>ij}?=i) ~ A - If a number £ £ {1, . . . , \S\} and 
an arbitrary direction is assigned to each edge {i,j} £ £, the 
(oriented) incidence matrix B £ R nX l £ l is defined component- 
wise as Bki = 1 if node k is the sink node of edge £ and as 
Bk£ — — 1 if node k is the source node of edge £; all other 
elements are zero. For x £ R n , the vector B T x has com- 



ponents Xi — Xj for any oriented edge from j to z, that is, 
B T maps node variables Xi, Xj to incremental edge variables 
Xi — xj. If diag({aij}{ i; j} G £;) is the diagonal matrix of nonzero 
edge weights, then L = B dia,g({aij}{ijy e g)B T . For a vector 

x £ R n , the incremental norm ||x||^,oo = max{ijy e g used in 
the main article, can be expressed via the incidence matrix 
B as ||x||£,oo = ||^ T ^||oo. If the graph is connected, then 
Ker (B T ) — Ker (L) = span(l n ), all n — 1 remaining eigenval- 
ues of L are real and strictly positive, and the second-smallest 
eigenvalue \2(L) is called the algebraic connectivity. The or- 
thogonal vector spaces Ker (B) and Ker (B) ± = lm(B T ) are 
spanned by vectors associated to cycles and cut-sets in the 
graph , see for example [1, Section 4] or [2]. In the following, 
we refer to Ker (B) and Im (B T ) as the cycle space and the 
cut-set space, respectively. 

Laplacian inverses: Since the Laplacian matrix L is sin- 
gular, we will frequently use its Moore-Penrose pseudo in- 
verse . If U £ R nXn is an orthonormal matrix of eigen- 
vectors of L, the singular value decomposition of L is L = 
[7diag({0, A2, . . . , \n})U T ', and its Moore-Penrose pseudo in- 
verse is given by L f = C7diag({0, 1/A 2 , . . . , l/\ n })U T . We 
will frequently use the identity L • — • L — I n — ^l nX n, 
which follows directly from the singular value decomposition. 
We also define the effective resistance between nodes i and j 
by Rij = Lj. + L\. - 2L\.. We refer to [3] for further infor- 
mation on Laplacian inverses and on the resistance distance. 

Mathematical Models and Synchronization Notions 

In this section we introduce the mathematical model of cou- 
pled phase oscillators considered in this article, we present 
some synchronization notions, and give a detailed account of 
the literature on synchronization of coupled phase oscillators. 

General Coupled Oscillator Model. Consider a weighted, undi- 
rected, and connected graph G — (V,S,A) with n nodes 
V = {1, . . . , n}, partitioned node set V = Vi U V2 and edge 
set S induced by the adjacency matrix A £ R nXn . We assume 
that the graph G has no self-loops {i, i}, that is, an = for all 
i £ V. Associated to this graph, consider the following model 
of |Vi| > second-order Newtonian and IV2I > first-order 
kinematic phase oscillators 

Midi + DiOi = un - > _ a,ij sm($i -Oj), i€Vi, 

[1] 

DiOi = Ui — ._ l cnj sin(6i — Oj) , i £ V2, 

where 6i £ S 1 and 6i £ R 1 are the phase and frequency of os- 
cillator i £ V, uji £ R 1 and Di > are the natural frequency 
and damping coefficient of oscillator i £ V, and Mi > is iner- 
tial constant of oscillator i £ Vi . The coupled oscillator model 
[1] evolves on T n x R |Vl1 , and features an important symme- 
try, namely the rotational invariance of the angular variable 6. 
The interesting dynamics of the coupled oscillator model [1] 
arises from a competition between each oscillator's tendency 
to align with its natural frequency uji and the synchronization- 
enforcing coupling aij sin(0i — Oj) with its neighbors. 

As discussed in the main article, the coupled oscillator 
model [1] unifies various models proposed in the literature. 
For example, for the parameters Vi = and Di — 1 for all 
i £ V2, it reduces to the celebrated Kuramoto model [4, 5] 

Oi = uji - } sin(0i - Oj) , i £ {1, . . . , n} . [2] 

We refer to the review articles [6, 7, 8, 9, 10] for various the- 
oretic results on the Kuramoto model [2] and further syn- 
chronization applications in natural sciences, technology, and 



7 



social networks. Here, we present a detailed modeling of the 
spring oscillator network used as a mechanical analog in the 
main article, and we present a few power network models, 
which can be described by the coupled oscillator model [ 1 ] . 

Mechanical Spring Network. Consider the spring network il- 
lustrated in Figure 6 consisting of a group of n particles con- 
strained to rotate around a circle with unit radius. For sim- 
plicity, we assume that the particles are allowed to move freely 
on the circle and exchange their order without collisions. 

Each particle is characterized by its phase angle ft £ S 1 
and frequency ft GR, and its inert ial and damping coefficients 
are Mi > and Di > 0. The external forces and torques 
acting on each particle are (i) a viscous damping force DiOi 
opposing the direction of motion, (ii) a non-conservative force 
uji £ R along the direction of motion depicting a preferred 
natural rotation frequency, and (iii) an elastic restoring torque 
between interacting particles i and j coupled by an ideal elas- 
tic spring with stiffness aij > and zero rest length. The 
topology of the spring network is described by the weighted, 
undirected, and connected graph G — (V, S,A). 

To compute the elastic torque between the particles, we 
parametrize the position of each particle i by the unit vector 
Pi = [cos (ft) , sin(ft)] T £ S 1 C R 2 . The elastic Hookean en- 
ergy stored in the springs is the function E : T n — » R given 
up to an additive constant by 

m = ^ {iJ}e£ a f\\ Pi - P s 

— dij (l — cos (ft) cos (ft ) — sin (ft) sin (ft)) 

where we employed the trigonometric identity cos (a — /3) = 
cos a cos /3 + sin a sin /3 in the last equality. Hence, we obtain 
the restoring torque acting on particle i as 

Ti{0) = ~Wi m = ~ ^J=i mj sin( ^ " 0j) ■ 

Therefore, the network of spring-interconnected particles de- 
picted in Figure 6 obeys the dynamics 

MiOi+DiOi =Ui-^ aijsm(0i-0j), ie{l,...,n}. [3] 

In conclusion, the spring network in Figure 6 is a mechanical 
analog of the coupled oscillator model [1] with V2 = 0. 



^^^^^^^^^^^ 

Fig. 6. Mechanical analog of the coupled oscillator model [1]. 

Power Network Model. The coupled oscillator model [1] in- 
cludes also a variety of power network models. We briefly 
present different power network models compatible with the 
coupled oscillator model [1] and refer to [11, Chapter 7] for a 
detailed derivation from a higher order first principle model. 

Consider a connected power network with generators Vi 
and load buses V2. The network is described by the symmet- 



ric nodal admittance matrix Y £ C nXn (augmented with the 
generator transient reactances) . If the network is lossless and 
the voltage levels \Vi\ at all nodes i £ Vi U V2 are constant, 
then the maximum real power transfer between any two nodes 
i, j £ Vi U V2 is aij = \Vi \ • \ Vj\ • ^s(Yij), where $s(Yij) denotes 
the susceptance of the transmission line £ S. With this 

notation the swing dynamics of generator i are given by 

•r — 

MiOi + Di6i = P m ,i - 2^ dij sin(ft - ft ) , ieVi, [4] 

where ft £ S 1 and ft £ R 1 are the generator rotor angle and 
frequency, ft £ S 1 for j £ V2 are the voltage phase angles at 
the load buses, and P m ,i > 0, Mi > 0, and Di > are the 
mechanical power input from the prime mover, the generator 
inertia constant, and the damping coefficient. 

For the load buses V2 , we consider the following three load 
models illustrated in Figure 7. 

1) PV buses with frequency -dependent loads: All load 
buses are PV buses, that is, the active power demand and 
the voltage magnitude \ Vi\ are specified for each bus. The real 
power drawn by load i consists of a constant term P\^ > and 
a frequency dependent term D^ft with Di > 0, as illustrated 
in Figure 7(a). The resulting real power balance equation is 

v — > n 

Di6i + P hi = -y aij sin(ft - ft) , i £ V 2 . [5] 
A — 'j=i 

The dynamics [4]-[5] are known as structure-preserving 
power network model [12], and equal the coupled oscillator 
model [1] for uji — P m ,i, i G Vi, and uji — —P\,i, i £ V2. 

2) PV buses with constant power loads: All load buses 
are PV buses, each load features a constant real power de- 
mand P\,i > 0, and the load damping in [5] is neglected, 
that is, Di = in equation [5]. The corresponding circuit- 
theoretic model is shown in Figure 7(b). If the angular dis- 
tances |ft(t) —0j(t)\ < 7r/2 are bounded for each transmission 
line {i,j} £ S (this condition will be precisely established in 
the next section), then the resulting differential-algebraic sys- 
tem has the same local stability properties as the dynamics 
[4]- [5], see [13]. Hence, all of our results apply locally also 
to the structure-preserving power network model [4]-[5] with 
zero load damping Di — for i £ V2. 

3) Constant current and constant admittance loads: If 
each load i £ V2 is modeled as a constant current demand 
Ii and an (inductive) admittance Y{ } shunt to ground as illus- 
trated in Figure 7(c), then the linear current-balance equa- 
tions are / = YV, where / £ C n and V £ C n are the vectors 
of nodal current injections and voltages. After elimination of 
the bus variables Vi, i £ V2, through Kron reduction [3], the 
resulting dynamics assume the form [3] known as the (loss- 
less) network-reduced power system model [14, 15]. We refer to 
[11, 3] for a detailed derivation of the network-reduced model. 



(a) (b) (c) 




Fig. 7. Equivalent circuits of the frequency-dependent load model (a), the constant 
power load model (b), and the constant current and admittance load model (c). 

To conclude this paragraph on power network modeling, 
we remark that a first-principle modeling of a DC power 
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source connected to an AC grid via a droop-controlled inverter 
results also in equation [5]; see [16] for further details. 

Synchronization Notions. The following subsets of the n-torus 
T n are essential for the synchronization problem: For 7 G 
[0, 7r/2[, let Ag(j) C T n be the closed set of angle arrays 
(0i,..., n ) with the property \0i — 0j\ < 7 for {i,j} G S. 
Also, let Ag(j) be the interior of Ag(j)- 

Definition 1. A solution (0,0) : R> (T n ,R |Vl1 ) to the 
coupled oscillator model [1] is said to be synchronized if 
0(0) G Ag(j) and there exists u sy nc G M 1 suc/i t/iat 0(£) = 
0(0) + cjsynclnt (mod 27r) and 0(t) = cj sy ncl|Vi| for all t > 0. 
In other words, here, synchronized trajectories have the prop- 
erties of frequency synchronization and phase cohesiveness, 
that is, all oscillators rotate with the same synchronization 
frequency cj sy nc and all their phases belong to the set Ag(j)- 
For a power network model [4]- [5], the notion of phase co- 
hesiveness is equivalent to bounded flows \ai 3 ; sin(0i — 6j)\ < 
dij sin(7) for all transmission lines {i, j} G S. 

For the coupled oscillator model [ 1 ] , the explicit synchro- 
nization frequency is given by cj sync = J2i=i ^/ J2?=i D ^ 
see [9] for a detailed derivation. By transforming to a ro- 
tating frame with frequency cj syn c and by replacing Ui by 
uji — DiCJsync, we obtain cj sy nc = (or equivalent ly uj G l n ) 
corresponding to balanced power injections ^2 ieVl Pm,i = 
X^ev 2 m P ower network applications. Hence, without loss 
of generality, we assume that uj G 1^ such that cj sy nc = 0. 

Given a point r G S 1 and an angle s G [0, 27r], let 
rot s (r) G S 1 be the rotation of r counterclockwise by the angle 
s. For (n, . . . , r n ) G T n , define the equivalence class 

[(n,...,r n )] = {(rot s (n ),..., rot s (r n ) G T n | s G [0,2tt]}. 

Clearly, if (n, . . . ,r n ) G A G (7) 5 then [(n, . . . ,r n )] C A G (7)- 
Definition 2. Given G Ag(t) /or some 7 G [0, 7r/2[, £/ie set 
([0], 0| Vl |) C T n x R' Vl ' is a synchronization manifold of the 
coupled oscillator model [ 1 ] . 

Note that a synchronized solution takes value in a synchro- 
nization manifold due to rotational symmetry. For two first- 
order oscillators [2] the state space T 2 , the set A g (tt/2), as 
well as the synchronization manifold [0*] associated to an an- 
gle array 0* = (0* , 0J) G T 2 are illustrated in Figure 8. 




Fig. 8. Illustration of the state space T 2 , the set Ag{tt/2), the synchronization 
manifold [0*] associated to a point #* = (0*,0J) G Ag(tt/2), the tangent 
space at 0* , and the translation vector I2. 

Existing Synchronization Conditions. The coupled oscillator 
dynamics [ 1 ] , and the Kuramoto dynamics [ 2 ] for that mat- 
ter, feature (i) the synchronizing effect of the coupling de- 
scribed by the weighted edges of the graph G(V, £, A) and (ii) 
the de-synchronizing effect of the dissimilar natural frequen- 
cies uj G 1^ at the nodes. Loosely speaking, synchronization 
occurs when the coupling dominates the dissimilarity. Various 



conditions are proposed in the power systems and synchro- 
nization literature to quantify this tradeoff between coupling 
and dissimilarity. The coupling is typically quantified by the 
algebraic connectivity A 2 (L) [17, 15, 18, 19, 20, 21] or the 
weighted nodal degree deg^ = YPj=i a H I 22 ' 3 > 23 > 15 > 24 1> 
and the dissimilarity is quantified by either absolute norms 
||cj|| p or incremental (relative) norms ||B T C(;||p, where typi- 
cally p G {2,oo}. Sometimes, these conditions can be evalu- 
ated only numerically since they are state-dependent [17, 22] 
or arise from a non-trivial linearization process, such as the 
Master stability function formalism [20, 21, 25]. In general, 
concise and accurate results are only known for specific topolo- 
gies such as complete graphs [9, 26] linear chains [27, 28] and 
complete bipartite graphs [29] with uniform weights. 

For arbitrary coupling topologies only sufficient condi- 
tions are known [17, 15, 18, 22] as well as numerical in- 
vestigations for random networks [30, 19, 31, 32]. To best 
of our knowledge, the sharpest and provably correct syn- 
chronization conditions for arbitrary topologies assume the 

form X 2 (L) > K ~ ^i| 2 ) ^ see [ 15 > Theorem 4.4]. 

For arbitrary undirected, connected, and weighted, graphs 
G(V,S,A), simulation studies indicate that the known suffi- 
cient conditions [17, 15, 18, 22] are conservative estimates on 
the threshold from incoherence to synchrony, and every review 
article on synchronization concludes with the open problem of 
finding sharp synchronization conditions [7, 9, 6, 20, 21, 33]. 

Mathematical Analysis of Synchronization 

This section presents our analysis of the synchronization prob- 
lem in the coupled oscillator model [ 1 ] . 

An Algebraic Approach to Synchronization. Here we present 
a novel analysis approach that reduces the synchronization 
problem to an equivalent algebraic problem that reveals the 
crucial role of cycles and cut-sets in the graph topology. In 
a first analysis step, we reduce the synchronization problem 
for the coupled oscillator model [1] to a simpler problem, 
namely stability of a first-order model. It turns out that exis- 
tence and local exponential stability of synchronized solutions 
of the coupled oscillator model [1] can be entirely described 
by means of the first-order Kuramoto model [ 2 ] . 

Lemma 1. (Synchronization equivalence) Consider the 
coupled oscillator model [1] and the Kuramoto model [2]. 
The following statements are equivalent for any 7 G [0,7r/2[ 
and any synchronization manifold ([^], 0|v x |) C Ag(7)xR' Vi '. 

(i) [0] is a locally exponentially stable synchronization mani- 
fold the Kuramoto model [2]; and 
(ii) ([^],0|y 1 |) is a locally exponentially stable synchronization 
manifold of the coupled oscillator model [ 1 ] . 

If the equivalent statements (i) and (ii) are true, then, locally 
near their respective synchronization manifolds, the coupled 
oscillator model [1] and the Kuramoto model [2] together 
with the frequency dynamics -f^ — —M~ 1 D6 are topologi- 
cal^ conjugate. 

Loosely speaking, the topological conjugacy result means 
that the trajectories of the two plots in Figure 9 can be con- 
tinuously deformed to match each other while preserving pa- 
rameterization of time. Lemma 1 is illustrated in Figure 9, 
and its proof can be found in [9, Theorems 5.1 and 5.3]. 
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Fig. 9. The left plot shows the phase space dynamics of a network of n = 4 
second-order oscillators [3] with V2 = and Kuramoto-type coupling dij = K/n 
for all distinct i,j G V\ = {1,...,4} and for K G K. The right plot 
shows the phase space dynamics corresponding to first-order Kuramoto oscillators 
[2] together with the frequency dynamics -^j = — M — 1 DO. The natural fre- 
quencies UJi and the coupling strength K are chosen such that o; S ync — and 
K = 1.1 • max^^j^ ...,4} \ UJ i~ UJ j\- From the same initial configuration #(0) 
(denoted by ■) both first and second-order oscillators converge exponentially to the 
same synchronized equilibria (denoted by as predicted by Lemma 1. 

By Lemma 1, the local synchronization problem for the 
coupled oscillator model [1] reduces to the synchronization 
problem for the first-order Kuramoto model [2]. Hence- 
forth, we restrict ourself to the Kuramoto model [2]. The 
following result is known in the synchronization literature 
[18, 15] as well as in power systems, where the saturation 
of a transmission line is corresponds to a singularity of the 
load flow Jacobian resulting in a saddle node bifurcation 
[34, 35, 13, 12, 36, 37, 22, 17, 38, 39, 40, 41]. 
Lemma 2. (Stable synchronization in Ag(tt/2)) Consider 
the Kuramoto model [2] with a connected graph G(V,£, A), 
and let 7 G [0,7r/2[. The following statements hold: 

1 ) Jacobian: The Jacobian of the Kuramoto model evaluated 
at G T n is given by 

J (9) = -Bdia,g({aij cos(0i - 6j)} {i j }e£ )B T ; 

2) Stability: If there exists an equilibrium point #* G Ag(7), 
then it belongs to a locally exponentially stable equilibrium 
manifold [0*] G Ag (7); and 

3) Uniqueness: This equilibrium manifold is unique in 
Ag(t)- 

Proof Since we have that -J^- (uJi — ^22=i aik sm (^ — 
0*0) = -Ylk=i a ^cos(0i-0 k ) and ^f- (uJi - J2k=i a ^ sin(^ - 

Ok)) = CLij cos(6i—6j), the negative Jacobian of the right-hand 
side of the Kuramoto model [ 2 ] equals the Laplacian matrix of 
the connected graph G(V,S, A) where dij = aij cos(0i — 0j). 
Equivalently, in compact notation the Jacobian is given by 
J (9) — — B di&g({aij cos(6i — 9j)}{i^y e s)B T . This completes 
the proof of statement 1). 

The Jacobian J (9) evaluated at an equilibrium point 
#* G Ag(7) is negative semidefinite with rank n — 1. Its 
nullspace is l n and arises from the rotational symmetry of the 
right-hand side of the Kuramoto model [2], see Figure 8 for an 
illustration. Consequently, the equilibrium point 0* G Ag(j) 
is locally (transversally) exponentially stable. Moreover, the 
corresponding equilibrium manifold [0*] G Ag(j) is locally ex- 
ponentially stable. This completes the proof of statement 2). 

The uniqueness statement 3) follows since the right-hand 
side of [2] is a one-to-one function for G Ag(tt/2), see [37, 
Corollary 1]. ■ 

By Lemma 2, the problem of finding a locally stable syn- 
chronization manifold reduces to that of finding a fixed point 



0* G Ag(7) for some 7 G [0, 7r/2[. The fixed-point equations 
of the Kuramoto model [2] read as 

En 
_^ dij sin(0i — 9j) , i G {1, . . . , n} . [6] 

In a compact notation the fixed-point equations [6] are 

u = 5diag ({aij}{ij }e e) sin(B T 0) . [7] 

The following conditions show that the natural frequencies 
uj have to be absolutely and incrementally bounded and the 
nodal degree has to be sufficiently large such that fixed points 
of [6] exist. 

Lemma 3. (Necessary synchronization conditions) Con- 
sider the Kuramoto model [2] with graph G(V,S,A) and 
U G l n . Let 7 G [0, tt/2[, and define the weighted nodal degree 
degi = X^j=i a ij f or eac h node i G {1, . . . , n}. The following 
statements hold: 

1) Absolute boundedness: If there exists a synchronized 
solution G Ag(7) ; then 

deg- sin(7) > \u)i\ for all i G {1, . . . ,n} . [8] 

2) Incremental boundedness: If there exists a synchro- 
nized solution G Ag(7) ; then 

(deg-+deg J )sin(7) > for all {ij} G S . [9] 

Proof The first condition arises since sin(ft — 6j) G 
[— sin(7), sin(7)] for G Ag(7), and the fixed-point equation 
[6] has no solution if condition [8] is not satisfied. 

Alternatively, since uj G In, a multiplication of the fixed 
point equation [7] by the vector (e^ — e 7 -) G 1^ , for {z, j} G £, 
or equivalently a subtraction of the zth and jth fixed-point 
equation [6], yields the following equation for all {i,j} G S: 

En 
(aik sin(6i - k ) - ajk sm(0j - 0k)) • [10] 
k — l 

Again, equation [10] has no solution in Ag(7) if condition 
[9] is not satisfied. ■ 
In the following we aim to find sufficient and sharp con- 
ditions under which the fixed-point equations [7] admit a so- 
lution 0* G Ag(7). We resort to a rather straightforward 
solution ansatz. By formally replacing each term sin(0i — 0j) 
in the fixed-point equations [7] by an auxiliary scalar variable 
ipij, the fixed-point equation [7] is equivalently written as 

uj = Bdiag ({aij}{ijy e e) ^, [H] 
^ = sin(B T 0) , [12] 

where ip G M'^' is a vector with elements ipij. We will refer 
to equations [11] as the auxiliary-fixed point equation, and 
characterize their properties in the following theorem. 
Theorem 1. (Properties of the fixed point equations) 

Consider the Kuramoto model [2] with graph G(V,S,A) and 
uj G In, its fixed-point equations [7], and the auxiliary fixed- 
point equations [11]. The following statements hold: 

1) Exact solution: Every solution of the auxiliary fixed- 
point equations [11] is of the form 

1p = B T L^UJ + t^hom , [13] 

where the homogeneous solution ^hom G satisfies 
diag ({aij}{ijy e e) Aom G Ker (B). 

2) Exact synchronization condition: Let 7 G [0,7r/2[. 
The following three statements are equivalent: 
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(i) There exists a solution 0* G Ag(t) to the fixed-point 
equation [7]; 

(ii) There exists a solution G Ag(j) to 

B T L^u + V^hom = sin(B T 6) . [14] 

for some VWi G diag ({l/a^} {i)j}e ^) ker(B); and 

fmj There exists a solution ip G M'^' to £/ie auxiliary fixed- 
point equation [11] of £/ie /orm [13] satisfying the 
norm constraint Halloo < sin(7) ana 7 £/ie c^/c/e constraint 
arcsin('0) G Im(5 T ). 
If the three equivalent statements (i), (ii), and (Hi) 
are true, then we have the identities B T Q* — B T — 
arcsin^). Additionally, [6*] G Ag(j) is a locally expo- 
nentially stable synchronization manifold. 

Proof Statement 1): Every solution ip G to the 

auxiliary fixed-point equations [11] is of the form i\) — 
i^hom + ^pt 5 where ^hom is the homogeneous solution and 
ippt is a particular solution. The homogeneous solution 
satisfies B diag ({(^ij}{i,j}es) ^hom = n . One can eas- 
ily verify that -0 pt = B T iJ uo is a particular solution 1 , 
since J Bdiag({a^}{ i;J } G £:)^pt = B diag({ai j } { ij }e s)B T L^uj = 

LL^UJ (l n - ^l nX n)w CU. 

Statement 2), equivalence ((i) <^4> (ii)) : If there exists a 
solution 0* of the fixed-point equations [7], then 0* can be 
equivalent ly obtained from equation [12] together with the 
solution [13] of the auxiliary equations [11]. These two equa- 
tions directly give equation [14]. 

Equivalence ((ii) <^ (hi)) : For #* G Ac (7), we have from 
equation [14] that Halloo < sin (7) and arcsin^) = B T 6*, 
that is, arcsin(V0 G lm(B T ). Conversely, if the norm con- 
straint 1 1 ip\\ 00 < sin (7) and the cycle constraint arcsin^) G 
Im (B T ) are met, then equation [14] is solvable in Ag(7), 
that is, there is 0* G Ag(j) such that arcsin(^) = B J 0* . 
The local exponential stability of the associated synchroniza- 
tion manifold [0*] follows then directly from Lemma 2. ■ 

The particular solution B T L^lo to the auxiliary fixed- 
point equations [11] lives in the cut-set space Ker(B) ± and 
the homogenous solution ^hom lives in the weighted cycle 
space ^hom G diag ({l/cLij}{ij}es) Ker (B). As a conse- 
quence, by statement (iii) of Theorem 1, for each cycle in the 
graph, we obtain one degree of freedom in choosing the ho- 
mogeneous solution ^hom as well as one nonlinear constraint 
c T arcsin^) = 0, where c G ker(B) is a signed path vector 
corresponding to the cycle. 

Remark 1. (Comments on necessity) The cycle space 
Ker (B) of the graph serves as a degree of freedom to find 
a minimum oo-norm solution ^* to equations [11] via 

mir VeRl £ l IHI°o subject to u = Bdiag {{aij} { ij }e s) ip~ 

* [15] 

By Theorem 1, such a minimum oo-norm solution ip* nec- 
essarily satisfies Halloo < sin(7) so that an equilibrium 
#* G Ag(7) exists. Hence, the condition Halloo < sin(7) 
is an optimal necessary synchronization condition. 

The optimization problem [15] - the minimum oo-norm 
solution to an under-determined and consistent system of lin- 
ear equations - is well studied in the context of kinemati- 
cally redundant manipulators. Its solution is known to be 
non-unique and contained in a disconnected solution space 
[42, 43]. Unfortunately, there is no "a priori" analytic formula 
to construct a minimum oo-norm solution, but the optimiza- 
tion problem is computationally tractable via its dual problem 
max^^R^ u T lo subject to || diag ({aij}{i,j}e£) B T u\\i — 1. □ 



Synchronization Assessment for Specific Networks. In this 
subsection we seek to establish that the condition 
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is sufficient for the existence of locally exponentially stable 
equilibria in Ag(tt/2). More general, for a given level of phase 
cohesiveness 7 G [0,7r/2[ we seek to establish that the condi- 
tion 
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is sufficient for the existence of locally exponentially stable 
equilibria in Ag(7). Since the right-hand side of [17] is a 
concave function of 7 G [0, 7r/2[ that achieves its supremum 
value at 7* = 7r/2, it follows that condition [17] implies [16]. 

In the main article, we provide a detailed interpretation 
of the synchronization conditions [16] and [17] from vari- 
ous practical perspectives. Before continuing our theoretical 
analysis, we provide two further abstract but insightful per- 
spectives on the conditions [16] and [17]. 

Remark 2. (Interpretation of the sync condition) 

Graph-theoretic interpretation: With regards to the exact and 
state-dependent norm and cycle conditions in statement (iii) 
of Theorem 1, the proposed condition [17] is simply a norm 
constraint on the network parameters in cut-set space Im (B T ) 
of the graph topology, and cycle components are discarded. 

Circuit-theoretic interpretation: In a circuit or power net- 
work, the variable uj G M n corresponds to nodal power injec- 
tions. Let x G M) £ \ satisfy Bx — cj, then x corresponds to 
equivalent power injections along lines {i,j} G S. 2 Condition 
[161 can then be rewritten as ||B t L^Bx|| < 1. The ma- 

L J II II 00 

trix B T L^B G R^ x l £ l has elements (e n - e? n ) T L\e k n - e n ) 
for {i,j},{k,£} G £, its diagonal elements are the effective 
resistances Rij, and its off-diagonal elements are the network 
distribution (sensitivity) factors [44, Appendix 11A]. Hence, 
from a circuit-theoretic perspective condition [16] restricts 
the pair-wise effective resistances and the routing of power 
through the network similar to the resistive synchronization 
conditions developed in [22, 3, 23] □ 

As it turns out, the exact state-dependent synchronization 
conditions in Theorem 1 can be easily evaluated for the spars- 
est (acyclic) and densest (homogeneous) topologies and for 
"worst-case" (cut-set inducing) and "best" (identical) natural 
frequencies. For all of these cases the scalar condition [17] is 
sharp. To quantify a "sharp" condition in the following the- 
orem, we distinguish between exact (necessary and sufficient) 
conditions and tight conditions, which are sufficient in general 
and become necessary over a set of parametric realizations. 

Theorem 2. (Sync condition for extremal network 
topologies and parameters) Consider the Kuramoto model 
[2] with connected graph G(V,S,A) and w G 1^. Consider 
the inequality condition [17] for 7 G [0, 7r/2[. 
The following statements hold: 



^Likewise, it can also be shown that (B diag({a^j }^ jy^g)) ' oj as well as 
diag({aij }^ 1 B^ oj are other possible particular solutions. All of 

these solutions differ only by addition of a homogenous solution. Each one can 
be interpreted as solution to a weighted least squares problem, see [42]. Fur- 
ther solutions can also be constructed in a graph-theoretic way by a spanning- 
tree decomposition, see [2]. Our specific choice -0pt = B T uj has the prop- 
erty that ip-pt £ Im(£? T ) lives in the cut-set space, and it is the most useful 
particular solution in order to proceed with our synchronization analysis. 
2 Notice that x is not uniquely determined if the circuit features loops. 
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(Gl) Exact synchronization condition for acyclic graphs: 

Assume that G(V,S,A) is acyclic. There exists an expo- 
nentially stable equilibrium 0* £ Ag(t) if and only if con- 
dition [17] holds. Moreover, in this case we have that 
B T 6* = arcsin(5 T L t o;) £ A G (7); 

(G2) Tight synchronization condition for homogeneous 
graphs: Assume that G(V,£,A) is a homogeneous graph, 
that is, there is K > such that aij = K for all distinct 
i,j £ {l,...,n}. Consider a compact interval ft C R, 
and Ze£ C l n 6e the set of all vectors with components 
fti £ H /or a// z £ {1, . . . , n}. For all uj £ ft there exists 
an exponentially stable equilibrium 0* £ Ag(t) if and only 
if condition [17] holds; 

(G3) Exact synchronization condition for cut-set induc- 
ing natural frequencies: Let Qi, Q2 £ and Ze£ 
fid" fre t/ie se£ of bipolar vectors with components fti £ 
{^1,^2} for i £ {1, . . . , n}. For all uj £ Lft there exists 
an exponentially stable equilibrium 0* £ Ag(t) if and only 
if condition [17] holds. Moreover, ft induces a cut-set: if 
IQ2 — fti\ = sin (7), £/ien /or cj = Lfi we obtain the stable 
equilibrium 6* £ Ag(7) satisfying B T 6* — aicsin(B T ft) , 
that is, for all {i,j} £ S, \9* — 0j \ = if fti = ftj and 
|6>* -6>*| =7 i/n< / n^; and 

fG^y) Asymptotic correctness: in £/ie Zzrm£ cj — >• n , 
£/iere exists an exponentially stable equilibrium 0* £ 
Ag(7) if and only if condition [17] holds. More- 
over, for each component i £ we have that 
\im^ 0n (B T 6*) i /( S Lrcsin(B T L^uj)) i = 1. 



Proof Statement (Gl): For an acyclic graph we have that 
Ker (B) — 0. According to Theorem 1, there exists an equilib- 
rium #* £ Ag(7) if and only if condition [17] is satisfied. In 
this case, we obtain B T Q* — arcsin(B T L^cj). This completes 
the proof of statement (Gl). 

Statement (G2): In the homogeneous case, we have that 
L = K(nl n -l n xn) and L 1 = ^(l n -^l reX n), see [3, Lemma 
3.13]. Thus, the inequality condition [17] can be equivalently 
rewritten as sin (7) > \\B T • uj II = ||B t cj|| . Accord- 

v// — II II 00 Kn \\ lloo 

ing to [9, Theorem 4.1], the Kuramoto model [2] with ho- 
mogenous coupling aij = K features an exponentially stable 
equilibrium #* £ Ag(7), 7 £ [0, 7r/2[, for all uj £ ft if and only 
if the condition K > ||B t cj|| In is satisfied. This concludes 

II lloo ' 

the proof of statement (G2). 

Statement (G3): For notational convenience, let c = Qi — 
H 2 . Then, for uj £ Lft, we have that B T L ] uj = B T L ] Lft = 
B T ft is a vector with components {— c, 0, +c}. Now consider 
the solution ip — B T L^uj — B T ft to the auxiliary fixed point 
equations [11], and notice that arcsin(^) = arcsin(B T fi) 
has components {— arcsin(c), 0, + arcsin(c)}. In particular, 
we have that arcsin(^) £ Im(B T ), and the exact synchro- 
nization condition from Theorem 1 is satisfied if and only if 
|| V II 00 = c < sin (7), which corresponds to condition [17]. 
The cut-set property follows since B T 6* = arcsin(^) has 
components {— arcsin(c), 0, + arcsin(c)} = {—7, 0, +7}. This 
concludes the proof of statement (G3). 

Statement (G4)' Since lim^o (arcsin(x)/x) = 1 for x £ 
R, we obtain lim^^o™ (ar csin(B T L^uj)i/(B T L^uj)) . = 1 for 
each component i £ {1,...,|£|}. Thus, the cycle constraint 
arcsin(^) £ Im (B T ) is asymptotically met with ip — B T L^uj. 
In this case, the solution of equation [14] is obtained as 
B T 6* — B t L^uj, and we have that #* £ Ag(7) if and only 
if the norm constraint [17] is satisfied. 3 This concludes the 
proof of statement (G4) and Theorem 2. ■ 



Theorem 1 shows that the solvability of the fixed-point 
equations [7] is inherently related to the cycle constraints. 
The following lemma establishes feasibility of a single cycle. 
Lemma 4. (Single cycle feasibility) Consider the Kuramoto 
model [2] with a cycle graph G(V,S,A) and uj £ 1^. With- 
out loss of generality, assume that the edges are labeled by 
{z, i + 1} (mod n) for i £ {1, . . . , n} and Ker (B) = span(l n ). 
Define x £ 1^ and y £ R^o uniquely by x = B T L^uj and 
yi = a i;( i + i) (mod n) > for i £ {1, . . . , n}. Let 7 £ [0, tt/2[. 
The following statements are equivalent: 

(i) There exists a stable equilibrium #* £ Ag(7); and 
(ii) The function f : [A m in, A ma x] —> R with domain boundaries 

\ — sin(7) — xa j \ • sin(7) — xa 

Amin — max ^ — l - and A max - mm — 4^ — l - 

iE{l,...,n} y% ie{l,...,n} Vl 

and defined by /(A) = 5^r=i arcsm (^i + A^) satisfies 
/(A min ) < < /(A max I • 

If both equivalent statements 1) and 2) are true, then B T 6* — 
arcsin(x + \*y), where A* £ [A m i n ,A ma x] satisfies /(A*) = 0. 

Proof According to Theorem 1, there exists a stable 
equilibrium 0* £ Ag(7) if and only if there exists a solu- 
tion ip = x + An, A £ R, to the auxiliary fixed-point equations 
[11] satisfying the norm constraint Halloo < sin(7) and the 
cycle constraint arcsin(^) £ lm(B T ). 

Equivalently, since Ker (B) = span(l n ), there is A £ R 
satisfying the norm constraint \\x + An||oo < sin(7) < 1 
and the cycle constraint 1 J arcsin(x + yX) — 0. Equiv- 
alently, the function /(A) = 1 J arcsin(x + yX) features a 
Zero A £ [Amin 5 Amax] (corresponding to the cycle constraint), 
where the constraints on A m in and A ma x guarantee the norm 
constraints xi + yiX max < sin(j) and Xi + ^A m in > — sin (7) 
for alH £ {1, . . . , n}. Equivalently, by the intermediate value 
theorem and due to continuity and (strict) monotonicity of 
the function /, we have that /(A m in) < < /(A max ). Finally, 
if A* £ [Amin, Amax] is found such that /(A*) = 0, then, by 
Theorem 1, B T Q* — arcsin('0) = arcsin(x + A*n). ■ 

Lemma 4 offers a checkable synchronization condition for 
cycles, which leads to the following theorem. 
Theorem 3. ( Sync conditions for cycle graphs) Con- 
sider the Kuramoto model [2] with a cycle graph G(V,£,A) 
and uj £ 1^ . Consider the inequality condition [17] for 
7 £ [0, 7r/2[. The following statements hold. 

(CI) Exact sync condition for symmetric natural fre- 
quencies: Assume that uj £ 1^ is such that B T L^uj is a 
symmetric vector 4 . There is an exponentially stable equi- 
librium 0* £ Ag(7) if and only if condition [17] holds. 
Moreover, in this case B T 6* =arcsin(,B T L t cj). 

(C2) Tight sync condition for low-dimensional cycles: 
Assume the network contains n £ {3, 4} oscillators. Con- 
sider a compact interval QcM, and let ft £ R n be the set 
of vectors with components fti £ ft for all i £ {1, . . . , n}. 
For all uj £ Lft there exists an exponentially stable equilib- 
rium 0* £ Ag(j) if and only if condition [17] holds. 

(C3) General cycles and network parameters: In general 
for n > 5 oscillators, condition [16] does not guarantee 
existence of an equilibrium 0* £ Ag(tt/2). As a sufficient 
condition, there exists an exponentially stable equilibrium 
6* € A ( 7 ), 7 € [0,n/2[, if 



3 Of course, the limit u> — > O n also implies that the resulting equilib- 
rium 0* G Ag(0) corresponds to phase synchronization 9i — 6j for all 
i, j G {l,...,n}. The converse statement 0* G Aq(0) uj — O n is 

also true and its proof can be found in [9, Theorem 5.5]. 

4 A vector x G 1^ is symmetric if its histogram is symmetric, that is, up to 
permutation of its elements, x is of the form x = [— c, +c] T for n even and 
some vector c G R n ^ 2 and x — [ — c, 0, +c] T for n odd and some c G IR^ 77 ' -1 )/ 2 . 
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\B T L*u)\ 



loo max{ij} e £ ciij + mm^ijy e£ en 



•sin( 7 ). [18] 



Proof To prove the statements of Theorem 3 and to show 
the existence of an equilibrium 0* £ Ag(7), we invoke the 
equivalent formulation via the function /(A) as constructed 
in Lemma 4. In particular, we seek to prove the statement: 



Let A n 



- sin(7) — Xi 



and A n 



= max< G{ i,..., n} y . 
sm(7)-a^ ^ f unc tion / : [A m i n ,A 



mm 5 '\maxj 



R defined by /(A) = J27=i arcsin(xi + Xyi) satisfies 
/(A m in) < < /(Amax) (equivalently there is A* £ 
[Amin, Amax] such that /(A*) = 0) if and only if the con- 
dition ||x||oo < sin (7) is satisfied. 

Statement (CI): For a symmetric vector x — B t L^uj, 
all odd moments about the (zero) mean vanish, that is, 
Er=i x * 2p+1 = for p G N . Since the Taylor series of the 
arcsin about zero features only odd powers, we have /(0) = 

£r=i arcsin(x,) = £?=i E^o 2^( P f)^2 P+ i) = °- State " 

ment 1) follows then immediately from Lemma 4. 

Statement (C2): By statement (CI), statement (C2) is 
true if B T L^u) is symmetric. Statement (C2), can then be 
proved in a combinatorial fashion by considering all devia- 
tions from symmetry arising for three or four oscillators. In 
order to continue recall that arcsin(x) is a super-additive func- 
tion for x £ [0, 1] and a sub-additive function for x £ [—1,0], 
that is, arcsin(x) + arcsin(^) < arcsin(x + y) for x, y > and 
x + y < 1, arcsin(x) + arcsin(y) > arcsin(x + y) for x,y < 
and x + 2/ > —1, and arcsin(x) + arcsin(^) = arcsin(x + y) for 
x = y — 0. We now consider each case n £ {3,4} separately. 

Proof of sufficiency for n — 3: Assume that ||x||oo < 
sin (7). Since the case /(A = 0) = 1 J arcsin(x) = for a 
symmetric vector x £ K 3 is already proved, we consider now 
the asymmetric case /(A = 0) = l^arcsin(x) > (the proof 
of the case 1„ arcsin(x) < is analogous). Necessarily, it 
follows that at least two elements of x are negative: if one 
element of x is zero, say x\ = 0, then we fall back into the 
symmetric case X2 = — X3; on the other hand, if only one 
element is negative, say x± < and X2,X3 > 0, then we ar- 
rive at a contradiction since /(A = 0) = J27=i arcsin(x^) = 
— arcsin(x2 + X3) + arcsin(x2) + arcsin(xs) < due to super- 
additivity and since x\ — —X2 — X3. Hence, without loss of 
generality, let x — [a-\-b, — a, — b] T where a, b > 0. By assump- 
tion ||x|| oo < sin(7). It follows that a + b < sin (7), a < sin(7), 



- sin(7) — Xi 



< 0. 



b < sin(7), and A m in — max iG {i r .. )n} ;/ 

Due to super-additivity, /(A = 0) = 1„ arcsin(x) = 
arcsin(a + b) — (arcsin(a) + arcsin(fr)) > 0. Now we evaluate 
/(A) at the lower end of its domain [A m in, A ma x] and obtain 

/(Amin) = arcsin(a + b + 2/1 A min ) + arcsin(-a + y 2 Xmin) 

+ arcsin(-6 + 2/3 Amin) • [19] 

By the definition of A m in, at least one summand on the 
right-hand side of [19] equals —7. Furthermore, notice that 
the second and the third summand are negative, and the 
first summand satisfies arcsin(a + b + 2/1 A m i n ) > —7. If 
arcsin(a + b + 2/1 A m in) = —7, then clearly /(A m in) < 0. In 
the other case, arcsin(a + b + 2/1 A m i n ) > —7, it follows that 

/(Amin) < arcsin(a + b + 2/iA mm ) - 7 
v v ' 

<o 

-|- max { arcsin ( — a -\- £/2A m in ),arcsin(-6 + 2/3 Amin)} <0. 

<o 



Since /(A m in) < < /(0) < /(A max ), it follows from Lemma 
4 that there exists a stable equilibrium #* £ Ag(7). The 
sufficiency is proved for n = 3. 

Proof of sufficiency for n — 4: Assume that ||x||oo < 
sin(7). Without loss of generality, let argmax^ ... ;4 }{|xi|} 
be a singleton (otherwise x is necessarily symmetric), and let 
x £ 1^; be such that /(A = 0) = l^arcsin(x) > (the proof 
of the case 1„ arcsin(x) < is analogous). Necessarily, it fol- 
lows that at least two elements of x are negative: if only one 
element of x is negative, say x\ < and X2, £3, X4 > 0, then we 
arrive at a contradiction since /(A = 0) = EILi arcsin(xi) = 
— arcsin(x2 + X3 + £4) + arcsin(x2) + arcsin(xs) + arcsin(x4) 
is zero only in the symmetric case (for example, X2 = X3 = 
< xa = —xi) and strictly negative otherwise (due to super- 
additivity). If exactly one element of x is positive (and three 
are non-postive) , say x = [a + 6 + c, —a, —6, — c] T for a, 6, c > 
and a + 6 + c — \\x\\oo < sin(7), then, an analogous reasoning 
to the case n — 3 leads to /(A m in) < 0. 

It remains to consider the case of two positive and two 
negative entries. Without loss of generality let xi > X2 > > 
X3 > X4, where x\ 7^ — X4 and X2 7^ — X3 (this is the sym- 
metric case), X^r=i Xi ~ ^5 an< ^ IWI 00 — sm (7) by assump- 



tion. It follows that A n 



iE{l,...,n} 



n (7)~ 



< 0. 



Since /(A = 0) = 1^ arcsin(x) > and l„x = 0, it fol- 
lows from super-additivity that ||x||oo = max{xi,X2}, and 
the set argmax{xi, X2} must be a singleton (otherwise we 
arrive again at a contradiction or at the symmetric case). 
Suppose that ||x||oo = max{xi,X2} = £1, then necessarily 
\x 2 \ < \x3\ < \xa\ < \xi \ < sin(7). It follows that A m in < 0. 

Again, we evaluate the sum /(A m in) = J2t=i arcsin (x^ + 
2/iA m in). Notice that the last two summands arcsin (X3 + 
2/3 Amin) and arcsin (xi + 2/4A m in) are negative (since > 
X3 > X4 and Amin < 0), and the first two summands sat- 
isfy min{arcsin(xi + y± A m in), arcsin (x 2 + 2/2Amin)} > -7- If 
min{arcsin(xi + y\ A m in), arcsin(x2 + 2/2A m in)} = —7, we have 

/(Amin) = arcsin(x 3 + 2/3A m in) + arcsin(x 4 + 2/4 Amin) 

v y 

<0 

+ (-7 + max{arcsin(xi + yi A m in), arcsin (x 2 +2/2A m in)}) < 0. 



In case that min {arcsin (a + y\ A m in), arcsin (b+ y 2 A m in)} > —7, 
we obtain min ie { 3) 4}{arcsin(xi + yiX m in)} — —7 and 

/(Amin) < arcsin(xi + 2/iA min ) + arcsin(x 2 + 2/2A m in) 7 
+ max { arcsin (pa + ^A m in)} • 

Since \x 2 \ < \x%\ < \x±\ < \x\\ < sin (7), it readily follows 
that arcsin(xi + 2/1 A m i n ) — 7 < and arcsin(x2 + 2/2A m in) + 
max lG { 3 4} {arcsin(xi + ^A m in)} < 0. We conclude that 
/(Amin) < 0. Since /(A min ) < < /(0) < /(A max ), it fol- 
lows from Lemma 4 that there exists a stable equilibrium 
0* £ Ag(j)> The sufficiency is proved for n — 4. 

Proof of necessity for n £ {3, 4}: We prove the necessity 
by contradiction. Consider a compact cube Q = [— c, +c]' £ ' C 
R' £ ', where c > satisfies c > sin (7). Assume that for every 
x £ l n , even those satisfying II 

X 00 > c, there exists A £ R 
such that the cycle constraint 1 J arcsin(x + Xy) = and 
the norm constraint ||x + A2/H 00 < sin (7) are simultaneously 
satisfied. For the sake of contradiction, consider now the sym- 
metric case, where x £ 1^ has components x% £ {— c, +c, 0}. 
As proved in statement (CI), A* = uniquely solves the cycle 
constraint equation = /(A* = 0) = Y^i=i arcsin(x; + X*y) = 
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Sr=i arcsin(zbc) for any value of c £ [0, 1]. However, the norm 
constraint ||x + A*?/||oo = ||#||oo < sin (7) can be satisfied only 
if Ik II 00 < sin (7) < c. We arrive at a contradiction since we 
assumed II > c > sin(7). 

We conclude that, if x = B T L^uj is bounded within a com- 
pact cube Q = [-c, +c] l<s| C with c < sin(7), the con- 
dition [17] is also necessary for synchronization of all con- 
sidered parametric realizations of B T L^uj within this com- 
pact cube Q. For the compact set ft = Q n £ M n , it follows 
that the image B T o Lfl — B T fl equals the compact cube 
Q — [— (max wG Q uj — min^n uj) , +(max we Q uj— min^n uj)] n . 
Hence, the condition [17] is necessary for synchronization of 
all considered parametric realizations of uj in the compact set 
Lfl. This concludes the proof of statement (C2). 

Statement (C3): To prove the first part of statement (C3) 
we construct an explicit counterexample. Consider a cycle of 
length n > 5 with unit- weighed edges a^+i = 1, and let 

^ = "•[1 + ^3 -2 1-^ n _ 4 ] T , 

where a £ [0, 1]. For a < 1, these parameters satisfy the nec- 
essary conditions [8] and [9]. For the given parameters, we 
obtain the no n- symmetric vector x — B T L^uj given by 

x = B T tfuj = a- [-1 -1 1 ^l (n _ 3) ] T . [20] 

Notice that ||x||oo = ot < 1, x is non- symmetric, and x is the 
minimum oo-norm vector ip = x + Al n for A £ R. 

In the following^ we will show that there exists no equi- 
librium in lim 7 ^ 7r / 2 Ag(j) — Ag(tt/2). Consider the function 
/(A) = arcsin(l^x + Al n ) whose domain is centered symmet- 
rically around zero, that is, A max = — A m i n = lim 7 ^ 7r / 2 (sin(7) — 
a) = 1 — a. Notice that the domain of / vanishes as a f 1. 
For n — ► 00 we have that lim n ^oo /(0) = — arcsin(a) + 
lim n _).oo(^ — 3) • arcsin(a/ (n — 3)) = — arcsin(a) + a. Hence, 
as n — ► 00 and a | 1, we obtain /(0) = — \ + 1 < 0. Due 
to continuity of / with respect to a, n, A, we conclude that 
for n > 5 sufficiently large and a < 1 sufficiently large, 
there is no A* such that /(A*) = 0. Hence, the condition 
WB 1 L bJ Woo < 1 does generally not guarantee ex- 
istence of 0* £ Ag(tt/2) d Ag(tt/2). A second numerical 
counterexample will be constructed in Example 1 below. 

A sufficient condition for the existence of an equilib- 
rium 0* £ Ag(j) is Xi + Amin^/i < < Xj -|- A m ax2/i 

for each z £ {l,...,n}, which is equivalent to condition 
[18]. Indeed if condition [18] holds, we obtain /(A m i n ) = 
Sr=i arcsm ( x * + Amin^/z) as a sum of nonpositive terms and 
/(A max ) = £r=i arcsin(x; + A ma x2/i) as a sum of nonnegative 
terms. Since l^x = and generally x / n (otherwise we 
fall back in the symmetric case), at least one Xi is strictly 
negative and at least one x% is strictly positive, and it follows 
that /(A m in) < < /(Amax). The statement (C3) follows then 
immediately from Lemma 4. This concludes the proof. ■ 
In the following, define a patched network {G(V, 8, A), uj} 
as a collection of subgraphs and natural frequencies uj £ 1^, 
where (i) each subgraph is connected, (ii) in each subgraph 
one of the conditions (G1),(G2),(G3),(G4), (CI), or (C2) is 
satisfied, (Hi) the subgraphs are connected to another through 
edges {i,j} £ 8 satisfying \\(e l \ £ \ — e J ^) T iJ uj\ \ 00 < sin (7), and 
(iv) the set of cycles in the overall graph G(V, 8, A) is equal to 
the union of the cycles of all subgraphs. Since a patched net- 
work satisfies the synchronization condition [17] as well the 
norm and cycle constraints, we can state the following result. 
Corollary 4. (Sync condition for a patched network) 
Consider the Kuramoto model [2] with a patched network 
{G(V, 8, A), uj}, and let 7 £_[0,7r/2[. There is an exponen- 
tially stable equilibrium 0* £ Ag(j) if condition [17] holds. 



Example 1. (Numerical cyclic counterexample and 
its intuition) In the proof of Theorem 3, we provided an ana- 
lytic counterexample which demonstrates that condition [17] 
is not sufficiently tight for synchronization in sufficiently large 
cyclic networks. Here, we provide an additional numerical 
counterexample. Consider a cycle family of length n — 5 + 3 -p, 
where p £ No is a nonnegative integer. Without loss of gener- 
ality, assume that the edges are labeled by {z, i + 1} (mod n) 
for i £ {1, . . . , n} such that Ker (B) = span(l n ). Assume that 
all edges are unit-weighed a^+i ( mo d n) — 1 for i £ {1, . . . , n}. 
Consider a £ [0, 1[, and let 

uj = a • [-1/2 2 p+ i 3/2 2p +i] T . 

For n = 5 (p = 1) the graph and the network parameters are 
illustrated in Figure 10. For the given network parameters, 
we obtain the non-symmetric vector B t L^uj given by 

B T ]Juj = a-[l -1 -l( n _ 2 )/3 1/2 • l 2 ( n -2)/3] T • 

Analogously to the example provided in the proof of Theorem 
3, ||B t L^cj||oo = ot and B T L^uj is the minimum oo-norm vec- 
tor B T L^uj + Al n for A £ R. In the limit a t 1, the necessary 
condition [8] is satisfied with equality. In Figure 10, for a t 1, 
we have that UJ2 — 2, and the necessary condition [8] reads 
as ai2 + «23 = |^2| = 2, and the corresponding equilibrium 
equation sin(0i — 62) + sin(#3 — 62) — 2 can only be satisfied 
if 0i - 2 = tt/2 and 3 - 2 = tt/2. Thus, with two fixed 
edge differences there is no more "wiggle room" to compen- 
sate for the effects of uji, z £ {1,3,4,5}. As a consequence, 
there is no equilibrium 0* £ Ag(tt/2) for a = 1 or equivalently 
Ili^L^lloo = 1. Due to continuity of the equations [6] with 
respect to a, we conclude that for a < 1 sufficiently large 
there is no equilibrium either. Numerical investigations show 
that this conclusion is true, especially for very large cycles. 
For the extreme case p — 10 7 , we obtain the critical threshold 
a w 0.9475 where 0* £ A G (7r/2) ceases to exist. □ 




Fig. 10. Cycle graph with n = 5 nodes and non-symmetric choice of UJ. 



Notice that both the counterexample used in the proof 
of Theorem 3 and the one in Example 1 are at the bound- 
ary of the admissible parameter space, where the necessary 
condition [8] is marginally satisfied. In the next section, we 
establish that such "degenerate" counterexamples do almost 
never occur for generic network topologies and parameters. 

To conclude this section, we remark that the main tech- 
nical difficulty in proving sufficiency of the condition [17] for 
arbitrary graphs is the compact state space T n and the non- 
monotone sinusoidal coupling among the oscillators. Indeed, if 
the state space was R n and if the oscillators were coupled via 
non-decreasing and odd functions, then the synchronization 
problem simplifies tremendously and the counterexamples in 
the proof of Theorem 3 and in Example 1 do not occur; see 
[45] for an elegant analysis based on optimization theory. 
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Statistical Synchronization Assessment 

After having established that the synchronization condition 
[17] is necessary and sufficient for particular network topolo- 
gies and parameters, we now validate both its correctness and 
its accuracy for arbitrary networks. 

Statistical Assessment of Correctness. Extensive simulation 
studies lead us to the conclusion that condition [17] is correct 
in general and guarantees the existence of a stable equilibrium 
#* £ Ag(7). In order to validate this hypothesis we invoke 
probability estimation through Monte Carlo techniques, see 
[46, Section 9] and [47, Section 3] for a comprehensive review. 

We consider the following nominal random networks 
{G(V, £, A), uj} parametrized by the number n > 2 of nodes, 
the width a > of the sampling region for each natural fre- 
quency uji and i £ {1, . . . , n}, and a connected random graph 
model RGMO) = G(V,S(p)) with node set V = {l,...,n} 
and edge set £ — £{p) induced by a coupling parameter p £ 
[0, 1]. In particular, given the four parameters (n, RGM,p, a), 
a nominal random network is constructed as follows: 

(i) Network topology: To construct the network topology, we 
consider three different one-parameter families of random 
graph models RGM(p) = G(V,£(p)), each parameterized 
by the number of nodes n > 2 and a coupling parame- 
ter p £ [0, 1]. Specifically, we consider (i) an Erdos-Renyi 
random graph model (RGM = ERG) with probability p of 
connecting two nodes, (ii) a random geometric graph model 
(RGM = RGG) with sampling region [0, l] 2 C R 2 , con- 
nectivity radius p, and (iii) a Watts- Strogatz small world 
network (RGM = SMN) [48] with initial coupling of each 
node to its two nearest neighbors and rewiring probability 
p. If, for a given n > 2 and p £ [0, 1], the realization of a 
random graph model is not connected, then this realization 
is discarded and new realization is constructed; 

(ii) Coupling weights: For a given random graph G(V,£(p)), 
for each edge {i,j} £ £(p), the coupling weight = a,ji > 
is sampled from a uniform distribution supported on the 
interval [0.5, 5]; 

(iii) Natural frequencies: For a given n > 2 and a > 0, the nat- 
ural frequencies uj £ 1^ are constructed in two steps. In 
a first step, n real numbers qi, i £ {1, . . . , n}, are sampled 
from a uniform distribution supported on [— a/2,+a/2], 
where a > 0. In a second step, by subtracting the average 
Z)r=i we define Vi=qi- Sr=i for z £ {1, ... , n} 
and obtain uj — (cji, . . . , uj n ) £ l^"; and 

(iv) Parametric realizations: We consider forty realizations 
of the parameter 4-tuple (n,RGM,p,ct) covering a wide 
range of network sizes n, coupling parameters p, and natu- 
ral frequencies cj, which are listed in the first column of 
Table 3. The choices of a in these forty cases is such 
that the resulting equilibrium angles 0* satisfy on average 
max {i)j}G£ |0* - 0*| ^ tt/3. 

For each of the forty parametric realizations in (iv), we gen- 
erate 30000 nominal models of uj £ 1^ and G(V,£,A) (con- 
ditioned on connectivity) as detailed in (i) - (iii) above, each 
satisfying ||B T L^a;||cx3 < 1. If a sample does not satisfy 
Hi^L^cjUoo < 1, it is discarded and a new sample is gen- 
erated. Hence, we obtain 1.2 • 10 6 nominal random networks 
{G(y, £, A), cj}, each with a connected graph G(V,£,A) and 
uj £ 1^ satisfying ||5 t L"''cl;|| 00 < sin(7) for some 7 < tt/2. 

For each case and each instance, we numerically solve 
equation [7] with accuracy 10 -6 and test the hypothesis 

H:\\b t L^uj\\ <sin( 7 ) => 3 6>* £ A g (t) 

II 1 1 00 



with an accuracy 10" 4 . The results are reported in Table 3 
together with the empirical probability that the hypothesis H 
is true for a set of parameters (n,RGM,p,a). Given a set of 
parameters (n,RGM,p,a) and 30000 samples, the empirical 
probability is calculated as 

number of samples satisfying {% is true) 
Prob ( „, RGM , p , a) = . 

Given an accuracy level e £ ]0, 1[ and a confidence level 
77 £ ]0, 1[, we ask for the number of samples N such that the 
true probability Prob( n R GM,p,a) ((H is true) equals the empir- 
ical probability Prob( n)RGM)P)C ,) with confidence level greater 
than 1 — 77 and accuracy at least e, that is, 

Prob (jProb( n?RG M,p,cO (His true) - Prob (n)RG M,p,c) | < ej 

> 1 -77. 

By the Chernoff bound [46, Equation (9.14)], the number of 
samples N for a given accuracy e and confidence 77 is given as 

^>^°4 [21] 

For e = 7] = 0.01, the Chernoff bound [21] is satisfied for 
N > 26492 samples. By invoking the Chernoff bound [21], 
our simulations studies establish the following statement: 

With 99% confidence level, there is at least 99% ac- 
curacy that the hypothesis 7i is true with probability 
99.97 % for a nominal network constructed as in (i) - 
(iv) above. 

In particular, for a nominal network with parameters 
(n,RGM,p,a) constructed as in (i) - (iv) above, with 
99% confidence level, there is at least 99% accuracy that 
the probability Prob( n;R GM,p,a) (ji is true) equals the 

empirical probability Prob( n)R GM,p,«) 5 as listed in Ta- 
ble 3, that is, 



Prob^|Prob( n)RG M,p,c)( / H is true) 

— Prob( n;R GM,p,o;) I < 0.01^ > 0.99. 

It can be seen in Table 3 that for large and dense networks the 
hypothesis % is always true, whereas for small and sparsely 
connected networks the hypothesis 7i can marginally fail with 
an error of order O(10 -4 ). Thus, for these cases a tighter con- 
dition of the form ||B t L^cj||oo <sin(7) — O(10~ 4 ) is required 
to establish the existence of 0* £ Ag(7). These results in- 
dicate that "degenerate" topologies and parameters (such as 
the large and isolated cycles used in the proof of Theorem 3 
and in Example 1) are more likely to occur in small networks. 

Statistical Assessment of Accuracy. As established in the 
previous subsection, the synchronization condition [17] is 
a scalar synchronization test with predictive power for al- 
most all network topologies and parameters. This remark- 
able fact is difficult to establish via statistical studies in the 
vast parameter space. Since we proved in statement (G4) 
of Theorem 2 that condition [17] is exact for sufficiently 
small pairwise phase cohesiveness \0i — 6j\ <C 1 (or equiv- 
alently, for sufficiently identical natural frequencies uji and 
sufficiently strong coupling) , we investigate the other extreme 
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max{ijy e g \0i — 6j \ — 7r/2. To test the corresponding synchro- 
nization condition [16] in a low- dimensional parameter space, 
we consider a complex network of Kuramoto oscillators 

Oi = Ui - K • V n aij sin(0i -Oj), i e {1, . . . ,n} , [22] 

where K > is the coupling gain among the oscillators and 
the coupling weights are assumed to be unit- weighted, that is, 

= dji = 1 for all {i,j} G £. If L is the unweighted Lapla- 
cian matrix, then condition [16] reads as K > identical — 
||^^||^,oo. Of course, the condition K > identical is only suf- 
ficient and synchronization may occur for a smaller value of 
K than identical- In order to test the accuracy of the condi- 
tion K > if critical, we numerically found the smallest value of 
K leading to synchrony for various network sizes, connected 
random graph models, and sample distributions of the natu- 
ral frequencies. Here we discuss in detail the construction of 
the random network topologies and parameters leading to the 
data displayed in Figure 3 of the main manuscript. 

We consider the following nominal random networks 
{G(V, £, A), uj} parametrized by the number of nodes n £ 
{10,20,40,160}, the sampling distribution SD for the natu- 
ral frequencies uj £ 1^, and a connected random graph model 
RGM(p) = G(V, £(p)) with node set V = {1, ... , n} and edge 
set £ = £(p) induced by a coupling parameter p £ [0,1]. In 
particular, given the four parameters (n, RGM,p, SD), a nom- 
inal random network is constructed as follows: 

(i) Network topology and weights: To construct the network 
topology, we consider three different one-parameter fami- 
lies of random graph models RGM(p) = G(V,£(p)), each 
parameterized by the number of nodes n and a coupling 
parameter p £ [0,1]. Specifically, we consider (i) an Erdos- 
Renyi random graph model (RGM = ERG) with proba- 
bility p of connecting two nodes, (ii) a random geomet- 
ric graph model (RGM = RGG) with sampling region 
[0,1] 2 C M 2 , connectivity radius p, and (hi) a Watts- 
Strogatz small world network (RGM = SMN) [48] with 
initial coupling of each node to its two nearest neighbors 
and rewiring probability p. If, for a given n and p £ [0, 1], 
the realization of a random graph model is not connected, 
then this realization is discarded and new realization is 
constructed. All nonzero coupling weights are set to one, 
that is, = dji = 1 for {i,j} £ £; 

(ii) Natural frequencies: For a given network size n and sam- 
pling distribution SD, the natural frequencies uj £ 1^ are 
constructed in three steps. In a first step, the sampling 
distribution of the natural frequencies is chosen. For clas- 
sic Kuramoto oscillators with uniform coupling a%j — Kin 
for distinct i,j £ {l,...,n}, we know that the two ex- 
treme sampling distributions (with bounded support) are 
the bipolar discrete and the uniform distribution leading to 
the largest and smallest critical coupling, respectively [9]. 
Here we choose a uniform distribution (SD = uniform) 
supported on [— 1,+1] or a bipolar discrete distribution 
(SD = bipolar) supported on { — 1, +1}. In a second step, 
n real numbers qi, i £ {1, . . . , n}, are sampled from the dis- 
tribution SD. In a third step, by subtracting the average 
Sr=i Qi/ nwe define uji = qi - J27=i qi/n for i e {1, ... ,n} 
and obtain uj — (cji, . . . , uj n ) £ 1^; and 

(Hi) Parametric realizations: We consider 600 realizations of 
parameter 4-tuple (n, RGM,p, SD) covering a wide range 
of network sizes n, coupling parameters p, and natural fre- 
quencies uj. All 600 realizations are shown in Figure 3 in 
the main manuscript. 

For each of the 600 parametric realizations in (Hi), we gen- 
erate 100 nominal models of uj £ 1^ and G(V,£,A) (condi- 



tioned on connectivity) as detailed in (i) - (ii) above. Hence, 
we obtain 60000 nominal random networks {G(V, £, A), cj}, 
each with a connected graph G(V,£,A) and natural frequen- 
cies uj £ 1^. For each sample network, we consider the com- 
plex Kuramoto model [22] and numerically find the small- 
est value of K leading to synchrony with cohesive phases 
satisfying max^j}^ \0i — 0j\ = tt/2. The critical value of 
K is found iteratively by integrating the Kuramoto dynam- 
ics [22] and decreasing K if the steady state 0* satisfies 
max{y} € £ \0* — 0j \ < 7r/2 and increasing K otherwise. We 
repeat this iteration until a steady state 0* is found satisfying 
max{^} e £ \0i — 0j\ — 7r/2 with accuracy 10 -3 . Our findings 
are reported in Figure 3 in the main manuscript, where each 
data point corresponds to the sample mean of 100 nominal 
models with the same parameter 4-tuple (n, RGM,p, SD). 



Synchronization Assessment for Power Networks 

We envision that our proposed condition [17] can be applied 
to quickly assess synchronization and robustness in power net- 
works under volatile operating conditions. Since real- world 
power networks are carefully engineered systems with partic- 
ular network topologies and parameters, they cannot be re- 
duced to the standard topological random graph models [49], 
and we do not extrapolate the statistical results from the pre- 
vious section to power grids. Rather, we consider ten widely- 
established and commonly studied IEEE power network test 
cases provided by [50, 51] to validate the correctness and the 
predictive power of our synchronization condition [17]. 

Statistical Synchronization Assessment for IEEE Systems. 

We validate the synchronization condition [17] in a smart 
power grid scenario subject to fluctuations in load and gener- 
ation and equipped with fast-ramping generation and control- 
lable demand. Here, we report the detailed simulation setup 
leading to the results shown in Table 1 of the main manuscript. 

The nominal simulation parameters for the ten IEEE test 
cases can be found in [50, 51]. Under nominal operating con- 
ditions, the power generation is optimized to meet the forecast 
demand, while obeying the AC power flow laws and respecting 
the thermal limits of each transmission line. Thermal limits 
constraints are precisely equivalent to phase cohesiveness re- 
quirements, that is, for each line {i,j}, the angular distance 
\0i — 0j\ needs to be bounded such that the corresponding 
power flow a%j sin(6^ — 0j) is bounded. Here, we found the 
optimal generator power injections through the standard op- 
timal power flow solver provided by MATPOWER [50]. 

In order to test the synchronization condition [17] in a 
volatile smart grid scenario, we make the following changes to 
the nominal IEEE test cases with optimal generation: 

(i) Fluctuating loads with stochastic power demand: We as- 
sume fluctuating demand and randomize 50% of all loads 
(selected independently with identical distribution) to de- 
viate from the forecasted loads with Gaussian statistics 
(with nominal power injection as mean and standard devi- 
ation 0.3 in per unit system); 

(ii) Renewables with stochastic power generation: We assume 
that the grid is penetrated by renewables with severely fluc- 
tuating power outputs, for example, wind or solar farms, 
and we randomize 33% of all generating units (selected in- 
dependently with identical distribution) to deviate from 
the nominally scheduled generation with Gaussian statis- 
tics (with nominal power injection as mean and standard 
deviation 0.3 in per unit system); and 
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Fig. 11. Time series of the RTS 96 dynamics for 141% load increase resulting in ||B T Ltc*;|| 00 = ||Z/tu;||£ j00 = 0.9995 < 1. Figure (a) depicts the angles 0i(t), 
Figure (b) shows the frequencies 9i(t), and Figure (c) depicts the angular distances \6i(t) — 0j(t)\ over transmission lines, where the red dashed curves correspond to 
the pairs {121, 325} and {223, 318}. The inserts show the power injections CJj, the phase space of the generator dynamics (0(t), 0(t)), and the stationary angles 6^. 




(iii) Fast-ramping generation and controllable loads: Following 
the paradigm of smart operation of smart grids [52], the 
fluctuations can be mitigated by fast-ramping generation, 
such as fast-response energy storage including batteries 
and flywheels, and controllable loads, such as large-scale 
server farms or fleets of plug-in hybrid electrical vehicles. 
Here, we assume that the grid is equipped with 10% fast- 
ramping generation (10% of all generators, selected inde- 
pendently with identical distribution) and 10% controllable 
loads (10% of all loads, selected independently with identi- 
cal distribution), and the power imbalance (caused by fluc- 
tuating demand and generation) is uniformly dispatched 
among these adjustable power sources. 



For each of the ten IEEE test cases with optimal genera- 
tor power injections, we construct 1000 random realizations of 
the scenario (i)-(iii) described above. For each realization, we 
numerically check for the existence of a solution #* £ Ag(t), 
7 £ [0, 7r/2[ to the AC power flow equations, the right-hand 



side of the power network dynamics [4]- [5], given by 

En 
a,ij sin(0i - 6j) , ieVi, 

JZl [23] 

P\,i = - ._ 1 a ij sin(0i - 6j) , i G V 2 . 

The solution to the AC power flow equations [23] is found 
via the AC power flow solver provided by MATPOWER [50]. 
Notice that, by Lemma 2, if such a solution 6* exists, then it 
is unique (up to rotational invariance) and also locally expo- 
nentially stable with respect to the power network dynamics 
[4]- [5]. Next, we compare the numerical solution 0* with 
the results predicted by our synchronization condition [17]. 
As discussed in Remark 2, a physical insightful and compu- 
tationally efficient way to evaluate condition [17] is to solve 
the sparse and linear DC power flow equations given by 

En 
aij(5i - Sj) , % G Vi , 

^ [24] 

P\,i = - y i aij(6i - Sj) , i e v 2 . 
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The solution 5* of the DC power flow equations [24] is de- 
fined uniquely up to the usual translational invariance. Given 
the solution 5* of the DC power flow equations [24], the left- 
hand side of our synchronization condition [17] evaluates to 
II^L^Hoo = IlL^H^oo = max{ i)J -} G £ \6* - 5*\. 

Finally, we compare our prediction with the numerical re- 
sults. If II^L^Hoo < sin( 7 ) for some 7 £ [0, 7r/2[, then 
condition [17] predicts that there exists a stable solution 

e A G (7), or alternatively e A G (arcsin(|| J B T L t CL;|| 00 )). To 
validate this hypothesis, we compare the numerical solution 
0* to the AC power flow equations [23] with our prediction 
0* G Ao(siTCsm(\\B T L^(jj\\oo)). Our findings and the detailed 
statistics are reported in Table 1 of the main manuscript. It 
can be observed that condition [17] predicts the correct phase 
cohesiveness \0* — Q*\ along all transmission lines G S 
with extremely high accuracy even for large-scale networks, 
such as the Polish power grid model featuring 2383 nodes. 

Simulation Data for IEEE Reliability Test System 96. The 

IEEE Reliability Test System 1996 (RTS 96) is a widely 
adopted and relatively large-scale power network test case, 
which has been designed as a benchmark model for power 
flow and stability studies. The RTS 96 is a multi-area model 
featuring 40 load buses and 33 generation buses, as illustrated 
in Figure 4 in the main manuscript. The network parameters 
and the dynamic generator parameters can be found in [51]. 

The quantities in the coupled oscillator model [ 1 ] cor- 
respond to the product of the voltage magnitudes at buses i 
and j as well the susceptance of the transmission line con- 
necting buses i and j. For a given set of power injections at 
the buses and branch parameters, the voltage magnitudes and 
initial phase angles were calculated using the optimal power 
flow solver provided by MAT POWER [50]. The quantities 
Wi, i £ V2, are the real power demands at loads, and uji, 

1 G Vi, are the real power injections at the generators, which 
were found through the optimal power flow solver provided by 
MATPOWER [50]. We made the following changes in order 
to adapt the detailed RTS 96 model to the classic structure- 
preserving power network model [4]- [5] describing the gen- 
erator rotor and voltage phase dynamics. First, we replaced 
the synchronous condenser in the original RTS 96 model [51] 
by a U50 hydro generator. Second, since the numerical val- 



ues of the damping coefficients Di are not contained in the 
original RTS 96 description [51], we chose the following val- 
ues to be found in [53]: for the generator damping, we chose 
the uniform damping coefficient Di = 1 in per unit system 
and for i £ Vi , and for the load frequency coefficient we chose 
Di = 0.1s for i G V2. Third and finally, we discarded an 
optional high voltage DC link for the branch {113, 316}. 

Bifurcation Scenario in the IEEE Reliability Test System 96. 

As shown in the main manuscript, an imbalanced power dis- 
patch in the RTS 96 network together with a tripped gener- 
ator (generator 323) in the Southeastern (green) area results 
in a loss of synchrony since the maximal power transfer is 
limited due to thermal constraints. This loss of synchrony 
can be predicted by our synchronization condition [17] with 
extremely high accuracy. In the following, we show that a 
similar loss of synchrony occurs, even if the generator 323 is 
not disconnected and there are no thermal limit constraints 
on the transmission lines. In this case, the loss of synchrony is 
due to a saddle node bifurcation at an inter- area angle of 7r/2, 
which can be predicted accurately by condition [17] as well. 

For the following dynamic simulation we consider again an 
imbalanced power dispatch: the demand at each load in the 
Southeast (green) area is increased by a uniform amount and 
the resulting power imbalance is compensated by uniformly 
increasing the generation at each generator in the two West- 
ern (blue) areas. The imbalanced power dispatch essentially 
transforms the RTS 96 into a two-oscillator network, and we 
observe the classic loss of synchrony through a saddle-node 
bifurcation [9, 40] shown in Figures 11 and 12. In particular, 
the network is still synchronized for a load increase of 141% re- 
sulting in ||L^o;|| = 0.9995 < 1. If the loads are increased 

by an additional 10% resulting in IIL^cjIL = 1.0560 > 1, 

j 11 1 if., 00 

then synchronization is lost and the areas separate via the 
transmission lines {121,325} and {223,318}. Of course, in 
real-world power networks the transmission lines {121,325} 
and {223, 318} would be separated at some smaller inter-area 
angle 7* <C tt/2 due to thermal limits. For instance, the trans- 
mission line {121, 325} is separated at the angle 7* = 0.19787T 
corresponding to a 78% load increase, which can be predicted 
from condition [16] as 7* « arcsin(||L^a;||£ >00 ) with an accu- 
racy of 0.00237T. In summary, this transmission line scenario 
illustrates the accuracy of the proposed condition [16]. 
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Table 3. Results of the Monte Carlo simulations to test the hypothesis 7i. 



nominal random network 


failures of hypothesis 1~L: 


empirical probability: 


parametrized by (n,RGM,p,a) 


# (?{ is not true) 


Prob( n? RGM,;p,c0 


(10, ERG, 0.15, 6) 


104 


99.653 % 


(10, ERG, 0.3, 8) 


65 


99.783 % 


(10, ERG, 0.5, 14) 


15 


on oirr» 0/ 

99.950 % 


(10, ERG, 0.75, 25) 





100 % 


(20, ERG, 0.15, 10) 


on 

oO 


on 700 0/ 

99.733 % 


(^ZU, H/lHor, U.O, 10 ) 


D 


yy.yoo /o 


(90 ERG 5 24"i 

IZ/U, -L/J-lAJ, U.O, 


o 


100 % 


(90 F/RG 7^ 4^ 


n 
u 


-LUU /0 


(10 F/RG 1^ 1 1) 

yOU , -L/llAJ , U. -L O, J-oy 


22 


QQ Q97 0/ 


^OUj 1 J l v v J . U.O, ZU y 


n 
u 


1UU /o 


^OUj 1 j 1 1 vi . U.O, O I y 


n 
u 


100 % 
1UU /o 


(30 ERG 75 65 ^ 


o 


100 % 


(ftO F/RG 1 ^ 90^1 

^UU , -L/llAJ , U. J- O, ZiUy 


I 


QQ QQ7 0/ 
yy.yy / /o 


(60 F/RG 1 4fh 

^UU, -L/llAJ, U.O, ^±U J 


u 


100 % 

1UU /o 


(60 F/RG ^ 70) 

^UU, -L/IiAj, U.O, 1 U j 


o 


100 % 

-LUU /0 


(60, ERG, 0.75, 125) 





100 % 


(120, ERG, 0.15, 35) 





100 % 


(120, ERG, 0.3, 75) 





100 % 


(120, ERG, 0.5, 130) 





100 % 


(120, ERG, 0.75, 235) 





100 % 


(AO RGG ^ Ifh 


-LO 


QQ QRO % 

i7i7.i70U /o 


MO RGG 0^1^ 

^ ± U , 1 V v J u , U . O , J. O y 


1 ft 

-LO 


QQ Q40 % 

i7i7.i7 i +U /0 


(90 RGG 1 1 D^l 

^ZiU, llAJVJ, u.o, ±u y 


9*3 
^.o 


QQ Q94 0/ 


(90 RGG 0^1^ 

^ZiU, IiAjVj, U.O, -LO y 


-5 

o 


QQ QQ0 % 
yy.yyu /o 


(^o rgg n i 1 rh 

^ou, iiajvj, u.o, J-Uy 


31 


QQ RQ7 % 
yy.oj/ /o 


(30, RGG, 0.5, 15) 


1 


99.997 % 


(60, RGG, 0.3, 10) 


3 


99.990 % 


(60, RGG, 0.5, 15) 





100 % 


(120,RGG,0.3,10) 





100 % 


(120, RGG, 0.5, 15) 





100 % 


(10,SMN,0.1,10) 


2 


99.994 % 


(10,SMN,0.2,10) 





100 % 


(20,SMN,0.1,13) 





100 % 


(20,SMN,0.2,13) 





100 % 


(30,SMN,0.1,10) 





100 % 


(30,SMN,0.2,13) 





100 % 


(60,SMN,0.1,7) 





100 % 


(60,SMN,0.2,7) 





100 % 


(120,SMN,0.1,4) 





100 % 


(120,SMN,0.2,4) 





100 % 


over all 1.2 • 10 6 instances 


388 


99.968 % 



^Overall, 1.2 • 10 6 instances of {G(V, 8, A),oj} were constructed as described in (i) - (iv) above, each satisfying 
|| J B T L t c ( ;||oo < 1. For each instance, the fixed-point equation [7] was solved with accuracy 10 -6 , and the failures of the 
hypothesis H were reported within an accuracy of 10 ~ , that is, failures of order 10 ~ 5 were discarded. 
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